c     19/2/2000
c-----------------------------------------------------------------------------
c-----------------------------------------------------------------------------
c
c  MEBDF VERSION THAT USES YSMP (YALE SPARSE MATRIX PACKAGE)
c
c-----------------------------------------------------------------------------
c-----------------------------------------------------------------------------

      SUBROUTINE MEBDFSO(N,T0,HO,Y0,TOUT,TEND,MF,IDID,LOUT,
     +     LWORK, WORK,LIWORK,IWORK,MAXDER,ITOL,RTOL,ATOL,F,PDERV,
     +     NSP,NONZ,IPAR,RPAR,IERR)
c ****************************************************************************
c ****************************************************************************
c    Written by T.J. Abdulla and J.R. Cash,
c    Department of Mathematics,
c    Imperial College,
c    London SW7 2AZ
c    England
c    t.abdulla@ic.ac.uk   or  j.cash@ic.ac.uk
c    The authors would be pleased to receive any comments
c          good or bad!!
c   the authors are particulaly greatful to Torsten Hennig
c   from
c   Frauhoffer Institut fur Umweld-, Sicherheits- und Energietechnik,
c   UMSICHT,
c   Osterfelder Strasse 3
c   46047 Oberhausen
c   Germany
c      who provided us with a FORTRAN 90 version of this code.
c
c ****************************************************************************
c ****************************************************************************
c  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C     >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C     THIS SUBROUTINE IS FOR THE PURPOSE               *
C     OF SPLITTING UP THE WORK ARRAYS WORK AND IWORK   *
C     FOR USE INSIDE THE INTEGRATOR STIFF              *
C     <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C     .. SCALAR ARGUMENTS ..
      INTEGER IDID,LIWORK,LOUT,LWORK,MF,N,MAXDER,ITOL
C     ..
C     .. ARRAY ARGUMENTS ..
      DIMENSION  WORK(LWORK),Y0(N),RTOL(1),ATOL(1),
     &     IWORK(LIWORK),IPAR(*),RPAR(*)
C     ..
C     .. LOCAL SCALARS ..
      INTEGER I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13,
     +     JK1,JK2,JK3,JK4,JK5,JK6,JK7,JK8
C
      LOGICAL ORDER
C     COMMON BLOCKS
      COMMON /ORDERING/ORDER
C     ..
C     .. EXTERNAL SUBROUTINES ..
      EXTERNAL OVDRIV,F,PDERV
C     ..
C     .. SAVE STATEMENT ..
      SAVE  I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13,
     +     JK1,JK2,JK3,JK4,JK5,JK6,JK7,JK8
C      
      IF (IDID.EQ.1) THEN
         IF (N.LE.0) THEN
            WRITE (LOUT,9020) N
            IDID = -4
         ELSE
            I1 = N*12+ 3
            I2 = I1 + N*12
            I3 = I2 + N*2
            I4 = I3 + N
            I5 = I4 + N
            I6 = I5 + N
            I7 = I6 + N
            I8 = I7 + N
            I9 = I8 + N
            I10 = I9 + N
            I11 = I10 + NONZ
            I12 = I11 + NONZ
            I13 = I12 + NSP
C     
C.....INTEGER WORKSPACE
C     
            JK1 = N+15
            JK2 = JK1 + NONZ
            JK3 = JK2 + N
            JK4 = JK3 + N
            JK5 = JK4 + NONZ
            JK6 = JK5 + N
            JK7 = JK6 + N
            JK8 = JK7 + N
            UROUND=DLAMCH('Epsilon')
            WORK(1)=UROUND
            EPSJAC=SQRT(WORK(1))
C
            IF (LWORK.LT.(I13-1)) THEN
               IDID = -11
               WRITE (LOUT,9000) I13 - 1  
            ENDIF
            
            IF (LIWORK.LT.JK8-1) THEN
               IDID = -12
               WRITE (LOUT,9010) JK8 - 1
            ENDIF  
         ENDIF
         IF (IDID.LT.0) RETURN  
         ORDER = .FALSE.
      ENDIF
C
C    THE DIMENSION OF THE REAL WORKSPACE, WORK, HAS TO BE AT LEAST
C    (33*N + 2*NONZ + NSP + 3). WHILE THE DIMENSION OF THE INTEGER
c    WORKSPACE HAS TO BE AT LEAST (6*N + 2*NONZ + 15), WHERE NONZ IS 
C    THE NUMBER OF THE NONZERO ENTRIES IN THE JACOBIAN MATRIX (IF THIS 
C    IS UNKNOWN, USE AN ESTIMATE) AND NSP SHOULD BE LARGER THAN 
C    8*N+2+2*NONZ (WORKSPACE FOR THE SPARSE SOLVER) 
C
      CALL OVDRIV(N,T0,HO,Y0,TOUT,TEND,MF,IDID,LOUT,WORK(3),WORK(I1),
     +     WORK(I2),WORK(I3),WORK(I4),WORK(I5),WORK(I6),WORK(I7),
     +     WORK(I8),WORK(I9),WORK(I10),WORK(I11),WORK(I12),IWORK(13),
     +     IWORK(JK1),IWORK(JK2),IWORK(JK3),WORK(I12),MAXDER,ITOL,
     +     RTOL,ATOL,F,PDERV,IWORK(1),IWORK(2),IWORK(3),IWORK(4),
     +     IWORK(5),IWORK(6),IWORK(7),IWORK(8),IWORK(9),IWORK(10),
     +     IWORK(11),IWORK(12),UROUND,WORK(2),EPSJAC,NSP,NONZ,
     +     IWORK(JK4),IWORK(JK5),IWORK(JK6),IWORK(JK7),IPAR,RPAR,IERR)      

C     >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C     WORK() HOUSES THE FOLLOWING ARRAYS
C     Y(N,12)   , YHOLD(N,12) , YNHOLD(N,2) , YMAX(N)
C     ERRORS(N) , SAVE1(N) ,SAVE2(N),SCALE(N),ARH(N),PW(NONZ)
C     PWCOPY(NONZ) 
C     <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  
      RETURN
 9000 FORMAT (/,/,' ***** ERROR ***** INTEGRATION HALTED IN DRIVER',/,/,
     +       '   >>> REAL WORKSPACE IS INSUFFICIENT <<< ',/,
     +       '       WORKSPACE MUST BE AT LEAST ',I8,' ELEMENTS LONG')
 9010 FORMAT (/,/,' ***** ERROR ***** INTEGRATION HALTED IN DRIVER',/,/,
     +       '   >>> INTEGER WORKSPACE IS INSUFFICIENT <<< ',/,
     +       '       WORKSPACE MUST BE AT LEAST ',I8,' ELEMENTS LONG')
 9020 FORMAT (/,/,' ***** ERROR ***** INTEGRATION HALTED IN DRIVER',/,/,
     +       '   >>> ILLEGAL VALUE FOR NUMBER OF EQUATIONS <<< ',/,
     +       '                     WITH N = ',I8) 
      END 
C-----------------------------------------------------------------------------

      SUBROUTINE OVDRIV(N,T0,HO,Y0,TOUT,TEND,MF,IDID,LOUT,Y,YHOLD,
     +     YNHOLD,YMAX,ERRORS,SAVE1,SAVE2,SCALE,ARH,SAVF,PW,PWCOPY,
     +     RSP,IA,JA,P,IP,ISP,MAXDER,ITOL,RTOL,ATOL,F,PDERV,NQUSED,
     +     NSTEP,NFAIL,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET,MAXORD,MAXSTP,
     +     NNZ,UROUND,HUSED,EPSJAC,NSP,NONZ,IGP,JGP,INCL,JDONE,
     +     IPAR,RPAR,IERR) 
      
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C     THIS IS THE FEBRUARY 23rd 2000 VERSION OF OVDRIV, A PACKAGE FOR
C     THE SOLUTION OF THE INITIAL VALUE PROBLEM FOR SYSTEMS OF
C     ORDINARY DIFFERENTIAL EQUATIONS WITH SPARSE JACOBIANS
C     DY/DT = F(Y,T),    Y=(Y(1),Y(2),Y(3), . . . ,Y(N))
C     SUBROUTINE OVDRIV IS A DRIVER ROUTINE FOR THIS PACKAGE
C
C                    REFERENCES
C
C     1.  J. R. CASH, THE INTEGRATION OF STIFF INITIAL VALUE PROBLEMS 
C         IN O.D.E.S USING MODIFIED EXTENDED BACKWARD DIFFERENTIATION 
C         FORMULAE, COMP. AND MATHS. WITH APPLICS., 9, 645-657, (1983).
C     2.  J.R. CASH AND S. CONSIDINE, AN MEBDF CODE FOR STIFF
C         INITIAL VALUE PROBLEMS, ACM TRANS MATH SOFTWARE, 142-158,
C         (1992).
C     3.  J.R. CASH, STABLE RECURSIONS WITH APPLICATIONS TO THE
C         NUMERICAL SOLUTION OF STIFF SYSTEMS, ACADEMIC PRESS,(1979).
C     4.  A.C. HINDMARSH, ODEPACK, A SYSTEMISED COLLECTION OF ODE
C         SOLVERS, in SCIENTIFIC COMPUTING, R.S. STEPLEMAN et. al.
C         (eds) North-Holland, AMSTERDAM, pp55-64 , (1983).
C     5.  E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
C         EQUATIONS II, STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS,
C         SPRINGER 1996, page 267
C     6.  S.C. EISENSTAT, M.C. GURSKY, M.H. SCHULTZ, AND, A.H. SHERMAN,
C         YALE SPARSE MATRIX PACKAGE II, THE NONSYEMMETRIC CODES, RESEARCH 
C         REPORT NO. 114, DEPT. OF COMPUTER SCIENCES, YALE UNIVERSITY, 1977.
C-----------------------------------------------------------------------------
C-----------------------------------------------------------------------------
C     OVDRIV IS TO BE CALLED ONCE FOR EACH OUTPUT VALUE OF T, AND
C     IN TURN MAKES REPEATED CALLS TO THE CORE INTEGRATOR STIFF.
C
C     THE INPUT PARAMETERS ARE ..
C     N     =  THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS.
C     T0    =  THE INITIAL VALUE OF T, THE INDEPENDENT VARIABLE
C              (USED ONLY ON THE FIRST CALL)
C     HO    =  THE NEXT STEP SIZE IN T (USED FOR INPUT ONLY ON THE
C              FIRST CALL)
C     Y0    =  A VECTOR OF LENGTH N CONTAINING THE INITIAL VALUES OF Y
C              (USED FOR INPUT ONLY ON FIRST CALL)
C     TOUT  =  THE VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT.
C              INTEGRATION WILL NORMALLY GO SLIGHTLY BEYOND TOUT
C              AND THE PACKAGE WILL INTERPOLATE TO T = TOUT 
C     TEND  =  END OF THE RANGE OF INTEGRATION.
C     MF    =  THE METHOD FLAG.  AT PRESENT MF=25 OR MF = 26 ARE
C              ALLOWED. THESE ARE EXTENDED BACKWARD DIFFERENTIATION 
C              FORMULAE USING THE CHORD METHOD WITH ANALYTIC JACOBIAN
C              FOR MF = 25 AND NUMERICAL JACOBIAN FOR MF = 26.
C              WITH MF = 25 THE USER NEED  TO SPECIFY SUBROUTINE PDERV.
C     IDID  = THE INTEGER USED ON INPUT TO INDICATE THE TYPE OF CALL.
C              1   THIS IS THE FIRST CALL FOR THE PROBLEM.
C              0   THIS IS NOT THE FIRST CALL FOR THIS PROBLEM
C                  AND INTEGRATION IS TO CONTINUE.
C             -1   THIS IS NOT THE FIRST CALL FOR THE PROBLEM,
C                  AND THE USER HAS RESET N, RTOL, ATOL, AND/OR MF.
C              2   SAME AS 0 EXCEPT THAT TOUT HAS TO BE HIT 
C                  EXACTLY (NO INTERPOLATION IS DONE).
C                  ASSUMES TOUT .GE. THE CURRENT T.
C              3   SAME AS 0 EXCEPT CONTROL RETURNS TO CALLING
C                  PROGRAM AFTER ONE STEP. TOUT IS IGNORED, UNTIL THE 
C                  INTEGRATION REACHES TOUT OR BEYOND. IF IT PASSES TOUT
C                  THE PROGRAM INTERPOLATES THE SOLUTION VALUES AND
C                  RETURNS THE SOLUTION VALUE AT TOUT.
C              SINCE THE NORMAL OUTPUT VALUE OF IDID IS 0, 
C              IT NEED NOT BE RESET FOR NORMAL CONTINUATION.
C     LOUT   = THE LOGICAL OUTPUT CHANNEL FOR MESSAGE PASSING.
C     MAXDER=  THE MAXIMUM ORDER IS MAXDER + 1.
C              THE VALUE OF MAXDER CANNOT EXCEED 7.  THIS IS THE
C              VALUE RECOMMENDED UNLESS IT IS BELIEVED THAT THERE
C              ARE SEVERE STABILITY PROBLEMS IN WHICH CASE MAXDER=3
C              OR 4 SHOULD BE TRIED.
C     ITOL  =  AN INDICATOR OF THE TYPE OF ERROR CONTROL. SEE
C              DESCRIPTION BELOW UNDER ATOL. 
C     RTOL  =  A RELATIVE ERROR TOLERANCE PARAMETER. CAN BE EITHER A
C              SCALAR OR AN ARRAY OF LENGTH N.  SEE DESCRIPTION
C              BELOW UNDER ATOL.
C     ATOL  =  THE ABSOLUTE ERROR BOUND.
C              THE INPUT PARAMETERS ITOL, RTOL AND ATOL DETERMINE
C              THE ERROR CONTROL PERFORMED BY THE SOLVER.  THE 
C              SOLVER WILL CONTROL THE VECTOR e = (e(i)) OF ESTIMATED
C              LOCAL ERRORS IN y ACCORDING TO AN INEQUALITY OF THE FORM
C                  RMS-NORM OF (e(i)/ewt(i)) .LE. 1
C              THE ROOT MEAN SQUARE NORM IS
C                   RMS-NORM(V) = SQRT((SUM v(i)**2)/N).  HERE 
C                ewt = (ewt(i)) IS A VECTOR OF WEIGHTS WHICH MUST
C              ALWAYS BE POSITIVE, AND THE VALUES OF RTOL AND ATOL
C              SHOULD BE NON-NEGATIVE. IF ITOL = 1 THEN SINGLE STEP ERROR 
C              ESTIMATES DIVIDED BY YMAX(I) WILL BE KEPT LESS THAN 1
C              IN ROOT-MEAN-SQUARE NORM.  THE VECTOR YMAX OF WEIGHTS IS
C              COMPUTED IN OVDRIV. INITIALLY YMAX(I) IS SET AS
C              THE MAXIMUM OF 1 AND ABS(Y(I)).  THEREAFTER YMAX(I) IS
C              THE LARGEST VALUE OF ABS(Y(I)) SEEN SO FAR, OR THE
C              INITIAL VALUE YMAX(I) IF THAT IS LARGER.
C              IF ITOL = 1 THE USER NEEDS TO SET ATOL = RTOL =
C              THE PRECISION REQUIRED.  THEN
C                         ewt(i) = RTOL(1)*YMAX(i)
C              IF ITOL IS GREATER THAN 1 THEN
C                 ewt(i) = rtol(i)*abs(y(i)) + atol(i)
C              THE FOLLOWING TABLE GIVES THE TYPES (SCALAR/ARRAY)
C              OF RTOL AND ATOL, AND THE CORRESPONDING FORM OF ewt(i)
C                  ITOL   RTOL      ATOL       ewt(i)
C                   2    SCALAR    SCALAR   rtol*abs(y(i))   + atol
C                   3    SCALAR    ARRAY    rtol*abs(y(i))   + atol(i)
C                   4    ARRAY     SCALAR   rtol(i)*abs(y(i))+ atol
C                   5    ARRAY     ARRAY    rtol(i)*abs(y(i))+ atol(i)   
C              IF EITHER OF THESE PARAMETERS IS A SCALAR, IT NEED
C              NOT BE DIMENSIONED IN THE USER'S CALLING PROGRAM.
C
C     RPAR, IPAR  REAL AND INTEGER PARAMETERS (OR PARAMETER ARRAYS)
C              WHICH CAN BE USED FOR COMMUNICATION BETWEEN  THE USER'S
C              CALLING PROGRAM AND THE F, PDERV AND MAS SUBROUTINES.
C     IERR     IERR IS AN INTEGER FLAG WHICH IS ALWAYS EQUAL TO ZERO
C              ON INPUT.  SUBROUTINES F, PDERV AND MAS SHOULD ALTER
C              IERR ONLY IF ONE OF THEM ENCOUNTERS AN ILLEGAL OPERATION SUCH
C              AS THE SQUARE ROOT OF A NEGATIVE NUMBER OR EXPONENT
C              OVERFLOW. THE USER CAN THEN ALTER H AND CALL THE
C              SUBROUTINE AGAIN WITH IDID=-1 IF HE WISHES.
C
C
C     NONZ     THE NUMBER OF NONZERO ELEMENTS IN THE JACOBIAN MATRIX.
C              IF THIS UNKNOWN, USE AN ESTIMATE.
C
C     AFTER THE INITIAL CALL, IF A NORMAL RETURN OCCURED AND A NORMAL 
C     CONTINUATION IS DESIRED, SIMPLY RESET TOUT AND CALL AGAIN.
C     ALL OTHER PARAMETERS WILL BE READY FOR THE NEXT CALL. 
C     A CHANGE OF PARAMETERS WITH IDID = -1 CAN BE MADE AFTER
C     EITHER A SUCCESSFUL OR AN UNSUCCESSFUL RETURN.
C
C     THE OUTPUT PARAMETERS ARE..
C     T0    =  THE VALUE OF T WHICH RELATES TO THE CURRENT SOLUTION
C              POINT Y0()
C     HO    =  THE STEPSIZE H USED LAST, WHETHER SUCCESSFULLY OR NOT. 
C     Y0    =  THE COMPUTED VALUES OF Y AT T = TOUT
C     TOUT  =  UNCHANGED FROM ITS INPUT VALUE.
C     IDID  =  INTEGER USED ON OUTPUT TO INDICATE RESULTS, WITH
C              THE FOLLOWING VALUES AND MEANINGS..
C
C      0   INTEGRATION WAS COMPLETED TO TOUT OR BEYOND.
C
C     -1   THE INTEGRATION WAS HALTED AFTER FAILING TO PASS THE
C          ERROR TEST EVEN AFTER REDUCING H BY A FACTOR OF
C          1.E10 FROM ITS INITIAL VALUE.
C
C     -2   AFTER SOME INITIAL SUCCESS, THE INTEGRATION WAS
C          HALTED EITHER BY REPEATED ERROR TEST FAILURES OR BY
C          A TEST ON RTOL/ATOL.  TOO MUCH ACCURACY HAS BEEN REQUESTED.
C
C     -3   THE INTEGRATION WAS HALTED AFTER FAILING TO ACHIEVE
C          CORRECTOR CONVERGENCE EVEN AFTER REDUCING H BY A 
C          FACTOR OF 1.E10 FROM ITS INITIAL VALUE.
C
C     -4   IMMEDIATE HALT BECAUSE OF ILLEGAL VALUES OF INPUT
C          PARAMETERS.  SEE PRINTED MESSAGE.
C
C     -5   IDID WAS -1 ON INPUT, BUT THE DESIRED CHANGES OF
C          PARAMETERS WERE NOT IMPLEMENTED BECAUSE TOUT
C          WAS NOT BEYOND T.  INTERPOLATION AT T = TOUT WAS 
C          PERFORMED AS ON A NORMAL RETURN.  TO TRY AGAIN,
C          SIMPLY CALL AGAIN WITH IDID = -1 AND A NEW TOUT.
C
C
C     -6   MAXIMUM ALLOWABLE NUMBER OF INTEGRATION STEPS EXCEEDED.
C          TO CONTINUE THE USER SHOULD RESET IWORK(14).
C
C
C     -7   STEPSIZE IS TOO SMALL (LESS THAN SQRT(UROUND)/100)
C
C 
C     -11   INSUFFICIENT REAL WORKSPACE FOR THE INTEGRATION
C
C     -12   INSUFFICIENT INTEGER WORKSPACE FOR THE INTEGRATION
C
C     IN ADDITION TO OVDRIVE, THE FOLLOWING ROUTINES ARE PROVIDED
C     IN THE PACKAGE..
C
C     INTERP( - )   INTERPOLATES TO GET THE OUTPUT VALUES
C                   AT T=TOUT FROM THE DATA IN THE Y ARRAY. 
C     STIFF( - )    IS THE CORE INTEGRATOR ROUTINE.  IT PERFORMS A
C                   SINGLE STEP AND ASSOCIATED ERROR CONTROL.
C     COSET( - )    SETS COEFFICIENTS FOR BACKWARD DIFFERENTIATION
C                   SCHEMES FOR USE IN THE CORE INTEGRATOR. 
C     PSET( - )     COMPUTES AND PROCESSES THE JACOBIAN
C                   MATRIX J = DF/DY
C                   HAS BEEN CALLED FOR THE MATRIX A
C     SPARSDET( - ) CONSTRUCT THE SPARSITY STRUCUTRE OF THE JACOBIAN MATRIX
C                   AND COLUMN GROUPING FOR THE NUMERICAL JACOBIAN. 
C
C     AND THE FOLLOWING YALE SPARSE ROUTINES
C     OVDRV( - )    CONSTRUCT A REORDERING OF THE ROWS AND COLUMNS OF 
C                   THE ITERATION MATRIX BY MINIMUM DEGREE ALGORITHM
C     CDRV( - )     PERFORMS REORDERING, SYMBOLIC FACTORIZATION, NUMERICAL 
C                   FACTORIZATION OR LINEAR SYSTEM SOLUTION OPERATIONS 
C                   DEPENDING ON A PATH ARGUMENT.
C
C     THE FOLLOWING ROUTINES ARE TO BE SUPPLIED BY THE USER AND
C     SHOULD BE DECLARED AS EXTERNAL.                 
C
C     F(N,T,Y,YDOT,IPAR,RPAR,IERR)       
C                   COMPUTES THE FUNCTION YDOT = F(Y,T), THE
C                   RIGHT HAND SIDE OF THE O.D.E.
C                   HERE Y AND YDOT ARE VECTORS OF LENGTH N
C     
C     PDERV(N,T,Y,J,PD,IPAR,RPAR,IERR)   
C                   COMPUTES THE NONZERO ENTRIES IN THE Jth COLUMN OF THE 
C                   JACOBIAN MATRIX AND STORES IT IN PD, I.E., 
C                   PD(J) = DF(I)/DY(J) FOR ALL RELEVANT VALUES OF I.          
C                   PDERV IS CALLED ONLY IF MITER = 5. OTHERWISE A DUMMY
C                   ROUTINE CAN BE SUBSTITUTED.
C
C
C     THE DIMENSIONS OF YMAX,ERROR,SAVE1,SAVE2 AND THE FIRST DIMENSION
C     OF Y SHOULD ALL BE AT LEAST N.
C
C
C
C     UROUND   THIS IS THE UNIT ROUNDOFF AND HAS TO BE SET AS
C                        UROUND = DLAMCH('Epsilon')
C     EPSJAC   = SQRT(UROUND).
C
C     HUSED  (=WORK(2))    LAST STEPSIZE SUCCESSFULLY USED BY THE INTEGRATOR
C     NQUSED (=IWORK(1))   LAST ORDER SUCCESSFULLY USED
C     NSTEP  (=IWORK(2))   NUMBER OF SUCCESSFUL STEPS TAKEN SO FAR
C     NFAIL  (=IWORK(3))   NUMBER OF FAILED STEPS
C     NFE    (=IWORK(4))   NUMBER OF FUNCTION EVALUATIONS  SO FAR
C     NJE    (=IWORK(5))   NUMBER OF JACOBIAN EVALUATIONS  SO FAR
C     NDEC   (=IWORK(6))   NUMBER OF LU DECOMPOSITIONS  SO FAR
C     NBSOL  (=IWORK(7))   NUMBER OF 'BACKSOLVES'  SO FAR
C     NPSET  (=IWORK(8))   NUMBER OF TIMES A NEW COEFFICIENT MATRIX HAS BEEN
C               FORMED SO FAR
C     NCOSET (=IWORK(9))   NUMBER OF TIMES THE ORDER OF THE METHOD USED HAS
C               BEEN CHANGED SO FAR
C     MAXORD (=IWORK(10))  THE MAXIMUM ORDER USED SO FAR IN THE INTEGRATION
C     MAXSTP (=IWORK(11))  THE MAXIMUM ALLOWED NUMBER OF STEPS SET BY THE USER
C     NNZ    (=IWORK(12))  NUMBER OF NONZER ELEMENTS IN THE JACOBIAN MATRIX


C    IF IT IS ONLY REQUIRED TO CONTROL THE ACCURACY IN THE
C    DIFFERENTIAL VARIABLES THEN THE USER SHOULD FIND THE 
C    STRING 'AMMEND' AND MAKE CHANGES THERE
C
C     START OF PROGRAM PROPER 
C
C     .. SCALAR ARGUMENTS ..
      INTEGER IDID,LOUT,MF,N,MAXDER,ITOL,NZ
C     ..
C     .. ARRAY ARGUMENTS ..
      DIMENSION  ARH(N),ERRORS(N),PW(*),PWCOPY(*),SAVE1(N),
     +     SAVE2(N),Y(N,12),Y0(N),YHOLD(N,12),YMAX(N),YNHOLD(N,2),
     +     RTOL(1),ATOL(1),SCALE(N),RSP(NSP),IGP(NONZ),JGP(N),
     +     INCL(N),JDONE(N),IPAR(*),RPAR(*),SAVF(N) 
C     ..
C     .. LOCAL SCALARS ..
      INTEGER I,KGO,NHCUT,NSP,IA(N+1),JA(NONZ),IP(N),P(N),ISP(NSP)
C     ..
C     .. EXTERNAL SUBROUTINES ..
      EXTERNAL INTERP,STIFF,F,PDERV,SPARSDET
C     ..
C     .. INTRINSIC FUNCTIONS ..
      INTRINSIC DABS,DMAX1
C     ..
C     .. COMMON BLOCKS ..
      INTEGER IDOUB,ISAMP,JSTART,KFLAG,L,MAXORD,NBSOL,NCOSET,NDEC,
     +     NEWPAR,NFE,NJE,NPSET,NQUSED,NRENEW,NSTEP
C 
      SAVE T,H,HMIN,HMAX,KFLAG,JSTART
      LOGICAL ORDER
      COMMON/ORDERING/ORDER
C     ..
      IF (IDID.EQ.0) THEN
C        I.E. NORMAL CONTINUATION OF INTEGRATION 
         T0=T
         HMAX = DABS(TEND-T0)*10.0D+0
         IF ((T-TOUT)*H.GE.0.0D+0) THEN 
C           HAVE OVERSHOT THE OUTPUT POINT, SO INTERPOLATE
            CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0)
            IDID = KFLAG
            T0 = TOUT
            HO = H
            RETURN
         ENDIF
  
      ELSE IF (IDID.EQ.2) THEN
C        I.E. CONTINUING INTEGRATION BUT WISH TO HIT TOUT
         T0 = T
         HMAX = DABS(TEND-T0)*10.0D+0
         IF (((T+H)-TOUT)*H.GT.0.0D+0) THEN
C           WE HAVE ALREADY OVERSHOT THE OUTPUT POINT OR WE WILL
C           DO SO ON THE NEXT STEP
            IF (((T-TOUT)*H.GE.0.0D+0) .OR. (DABS(T-TOUT).LE.
     +           100.0D+0*UROUND*HMAX)) THEN
C              HAVE OVERSHOT THE OUTPUT POINT, SO INTERPOLATE
               CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0)
               T0 = TOUT
               HO = H
               IDID = KFLAG
               RETURN
  
            ELSE
C              WILL PASS TOUT ON NEXT STEP WITH CURRENT STEPSIZE
C              SO REDUCE STEPSIZE TO HIT TOUT 'EXACTLY'
               H = (TOUT-T)* (1.0D+0-4.0D+0*UROUND)
               JSTART = -1
            ENDIF
  
         ENDIF
  
      ELSE IF (IDID.EQ.-1) THEN
C        NOT FIRST CALL BUT PARAMETERS RESET
         H = HO
         IF(H.LT.EPSJAC/100.0D+0) THEN
            WRITE(LOUT,9160)
            IDID = -7
            RETURN
         ENDIF
         T0 = T
         IF ((T-TOUT)*H.GE.0.0D+0) THEN 
C           HAVE OVERSHOT TOUT
            WRITE (LOUT,9080) T,TOUT,H
            CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0)
            HO = H
            T0 = TOUT
            IDID = -5
            RETURN
  
         ELSE
            JSTART = -1
         ENDIF
  
      ELSE IF (IDID.EQ.3) THEN
         T0 = T
         IF ((T-TOUT)*H.GE.0.0D+0) THEN 
C           HAVE OVERSHOT,SO INTERPOLATE
            CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0)
            IDID = KFLAG
            T0 = TOUT
            HO = H
            RETURN
  
         ENDIF
         
      ELSE
C        IDID SHOULD BE 1 AND THIS IS THE FIRST CALL FOR THIS PROBLEM
C        CHECK THE ARGUMENTS THAT WERE PASSED FOR CORRECTNESS
         IF (IDID.NE.1) THEN 
C           VALUE OF IDID NOT ALLOWED
            WRITE (LOUT,9070) IDID 
            IDID = -4
         ENDIF
  
         NN=N
         IF(ITOL.LE.3) NN = 1
         DO 1 I=1,NN
            IF (RTOL(I).LT.0.0D+0) THEN
C           ILLEGAL VALUE FOR RELATIVE ERROR TOLERENCE
               WRITE (LOUT,9040)
               IDID = -4
            ENDIF
 1       CONTINUE

         NN=N
         IF(ITOL.EQ.1.OR.ITOL.EQ.2.OR.ITOL.EQ.4) NN=1
         DO 2 I=1,NN
            IF (ATOL(I).LT.0.0D+0) THEN
C           ILLEGAL ABSOLUTE ERROR TOLERANCE
               WRITE(LOUT,9045) 
               IDID=-4
            ENDIF
 2       CONTINUE
         IF(ITOL.EQ.1.AND.RTOL(1).EQ.0) THEN
C     ILLEGAL ERROR TOLERANCE
            WRITE(LOUT,9040)
            IDID = -4
         ENDIF
         IF(ITOL.NE.1) THEN
            VHOLD = 0.0D+0
            DO 3 I=1,N
               IF(ITOL.EQ.2) THEN
                  VHOLD = DMAX1(RTOL(1),ATOL(1))
               ELSE IF (ITOL.EQ.3) THEN
                  VHOLD = DMAX1(RTOL(1),ATOL(I))
               ELSE IF (ITOL.EQ.4) THEN
                  VHOLD = DMAX1(RTOL(I),ATOL(1))
               ELSE IF (ITOL.EQ.5) THEN
                  VHOLD = DMAX1(RTOL(I),ATOL(I))
               ENDIF
               IF(VHOLD.LE.0.0D+0) THEN
                  WRITE(LOUT,9040)
                  IDID = -4
               ENDIF
 3          CONTINUE
         ENDIF
         
         IF (N.LE.0) THEN
C     ILLEGAL VALUE FOR THE NUMBER OF EQUATIONS
            WRITE (LOUT,9050)
            IDID = -4
         ENDIF
         
         IF ((T0-TOUT)*HO.GE.0.0D+0) THEN
C           PARAMETERS FOR INTEGRATION ARE ILLEGAL
            WRITE (LOUT,9060)
            IDID = -4
         ENDIF
  
         IF ((MF.NE.25).AND.(MF.NE.26))THEN
C     ILLEGAL VALUE FOR METHOD FLAG
            WRITE (LOUT,9090) MF
            IDID = -4
         ENDIF
  
         IF(ITOL.LT.1.OR.ITOL.GT.5) THEN
C           ILLEGAL VALUE FOR ERROR CONTROL PARAMETER
            WRITE (LOUT,9110) 
            IDID=-4
         ENDIF

         IF(MAXDER.LT.1.OR.MAXDER.GT.7) THEN
C        ILLEGAL VALUE FOR MAXIMUM ORDER
            WRITE(LOUT,9120)
            IDID = -4
         ENDIF

         IF (IDID.NE.1) THEN 
            RETURN
         ELSE
C           THE INITIAL PARAMETERS ARE O.K. SO INITIALISE EVERYTHING
C           ELSE NECESSARY FOR THE INTEGRATION.
C           IF VALUES OF YMAX OTHER THAN THOSE SET BELOW ARE DESIRED, 
C           THEY SHOULD BE SET HERE. ALL YMAX(I) MUST BE POSITIVE. IF 
C           VALUES FOR HMIN OR HMAX, THE BOUNDS ON DABS(H), OTHER THAN
C           THOSE BELOW ARE DESIRED, THEY SHOULD BE SET BELOW.

            IF(ITOL.EQ.1) THEN
               DO 10 I = 1,N
                  YMAX(I) = DABS(Y0(I))
                  YMAX(I)=DMAX1(YMAX(I),1.0D0)
 10            CONTINUE
            ENDIF
            DO 15 I=1,N
               Y(I,1)=Y0(I)              
 15         CONTINUE
            T = T0
            H = HO
            HMIN = DABS(HO)
            HMAX = DABS(T0-TEND)*10.0D+0
            JSTART = 0
            NHCUT = 0
            
         ENDIF
      ENDIF

C
C     THE NEXT ROUTINE CONSTRUCT THE SPARSITY STRUCTURE OF THE ITERATION 
C     MATRIX (ROW, COLUMN) INDICES. THIS NEED TO BE DONE ONLY ONCE 
C     PER INTEGRATION
C
      IF (.NOT.  ORDER ) THEN
         DO 5 I=1,N
            AYI = DABS(Y(I,1))
            IF(ITOL.EQ.1) THEN
               SCALE(I) = YMAX(I)
            ELSE IF(ITOL.EQ.2) THEN
               SCALE(I) = RTOL(1)*AYI + ATOL(1)
            ELSE IF(ITOL.EQ.3) THEN
               SCALE(I) = RTOL(1)*AYI + ATOL(I)
            ELSE IF(ITOL.EQ.4) THEN
               SCALE(I) = RTOL(I)*AYI + ATOL(1)
            ELSE IF(ITOL.EQ.5) THEN
               SCALE(I) = RTOL(I)*AYI + ATOL(I)
            ENDIF
 5       CONTINUE
C
         CALL SPARSDET(N,T,Y0,PWCOPY,SCALE,IA,JA,SAVF,P,IP,NSP,
     +        ISP,RSP,NONZ,NGP,JGP,IGP,MF,INCL,JDONE,LOUT,F,PDERV,
     +        NNZ,IPAR,RPAR,IERR)                 
      ENDIF
      

C     <<<<<<<<<<<<<<<<<
C     <  TAKE A STEP  >
C     <<<<<<<<<<<<<<<<<
 20   IF ((T+H).EQ.T) THEN
         WRITE (LOUT,9000)
      ENDIF
C    
      CALL STIFF(H,HMAX,HMIN,JSTART,KFLAG,MF,T,TOUT,TEND,Y,N,
     +     YMAX,ERRORS,SAVE1,SAVE2,SCALE,PW,PWCOPY,YHOLD,YNHOLD,
     +     ARH,LOUT,MAXDER,ITOL,RTOL,ATOL,F,PDERV,NQUSED,NSTEP,
     +     NFAIL,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET,MAXORD,UROUND,
     +     EPSJAC,HUSED,RSP,IA,JA,P,IP,NSP,ISP,NONZ,SAVF,NGP,JGP,
     +     IGP,IPAR,RPAR,IERR)
C     
      IF (IERR .NE. 0)THEN
         WRITE(LOUT,9150)
         RETURN
      ENDIF
C     
      KGO = 1 - KFLAG
      IF (KGO.EQ.1) THEN
C        NORMAL RETURN FROM STIFF
         GO TO 30
         
      ELSE IF (KGO.EQ.2) THEN 
C        COULD NOT ACHIEVE REQUIRED PRECISION WITH HMIN
C        SO CHOP HMIN IF WE HAVEN'T DONE SO 10 TIMES
         GO TO 60
  
      ELSE IF (KGO.EQ.3) THEN 
C        ERROR REQUIREMENT SMALLER THAN CAN BE HANDLED FOR THIS PROBLEM 
         WRITE (LOUT,9010) T,H
         GO TO 70
  
      ELSE IF (KGO.EQ.4) THEN 
C        COULD NOT ACHIEVE CONVERGENCE WITH HMIN
         WRITE (LOUT,9030) T
         GO TO 60
      ENDIF
  
 30   CONTINUE
C ---------------------------------------------------------------------
C     NORMAL RETURN FROM THE INTEGRATOR.
C
C     THE WEIGHTS YMAX(I) ARE UPDATED IF ITOL=1.
C     IF DIFFERENT VALUES ARE DESIRED, THEY SHOULD BE SET HERE.
C
C     ANY OTHER TESTS OR CALCULATIONS THAT ARE REQUIRED AFTER EVERY
C     STEP SHOULD BE INSERTED HERE.
C
C     IF IDID = 3, Y0 IS SET TO THE CURRENT Y VALUES ON RETURN.
C     IF IDID = 2, H IS CONTROLLED TO HIT TOUT (WITHIN ROUNDOFF
C     ERROR), AND THEN THE CURRENT Y VALUES ARE PUT IN Y0 ON RETURN.
C     FOR ANY OTHER VALUE OF IDID, CONTROL RETURNS TO THE INTEGRATOR 
C     UNLESS TOUT HAS BEEN REACHED.  THEN INTERPOLATED VALUES OF Y ARE
C     COMPUTED AND STORED IN Y0 ON RETURN.
C     IF INTERPOLATION IS NOT DESIRED, THE CALL TO INTERP SHOULD BE
C     REMOVED AND CONTROL TRANSFERRED TO STATEMENT 500 INSTEAD OF 520.
C --------------------------------------------------------------------
      IF(NSTEP.GT.MAXSTP) THEN
         KGO=5
         KLAG=4
c   TOO MUCH WORK
         WRITE(LOUT,9130)
         IDID=-6
         GOTO 70
      ENDIF

      IF(ITOL.EQ.1) THEN
         D = 0.0D+0
         DO 40 I = 1,N 
            AYI = DABS(Y(I,1))
            YMAX(I)=DMAX1(YMAX(I),AYI)
 40      CONTINUE
      ENDIF
      IF (IDID.EQ.3.OR.IDID.EQ.1) GO TO 70
  
      IF (DABS(T-TOUT).LE.DABS(15.0D+0*UROUND*TOUT)) THEN
C        EFFECTIVELY WE HAVE HIT TOUT
         IDID = KFLAG
         T0 = TOUT
         DO 50 I = 1,N
            Y0(I) = Y(I,1)
 50      CONTINUE
         HO = H
         RETURN
         
      ENDIF
  
      IF (IDID.EQ.2) THEN
C        CONTINUING INTEGRATION BUT MUST HIT TOUT EXACTLY
         IF (((T+H)-TOUT)*H.GT.0.0D+0) THEN
C           WE HAVE ALREADY OVERSHOT THE OUTPUT POINT OR WE WILL DO
C           SO ON THE NEXT STEP
            IF (((T-TOUT)*H.GE.0.0D+0) .OR. (DABS(T-TOUT).LE.
     +           100.0D+0*UROUND*HMAX)) THEN
C              HAVE OVERSHOT, SO INTERPOLATE
               CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0)
               T0 = TOUT
               HO = H
               IDID = KFLAG
               RETURN
            ELSE
C              WILL PASS TOUT ON NEXT STEP WITH CURRENT STEPSIZE
C              SO REDUCE STEPSIZE TO HIT TOUT 'EXACTLY'
               H = (TOUT-T)* (1.0D+0-4.0D+0*UROUND)
               JSTART = -1
            ENDIF
         ENDIF
      ELSE IF ((T-TOUT)*H.GE.0.0D+0) THEN
C        HAVE OVERSHOT, SO INTERPOLATE
         CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0)
         IDID = KFLAG
         HO = H
         T0 = TOUT
         RETURN
      ENDIF
      
      GO TO 20
C ------------------------------------------------------------------- 
C     ON AN ERROR RETURN FROM THE INTEGRATOR, AN IMMEDIATE RETURN OCCURS
C     IF KFLAG = -2, AND RECOVERY ATTEMPTS ARE MADE OTHERWISE.
C     H AND HMIN ARE REDUCED BY A FACTOR OF .1 UP TO 10 TIMES
C     BEFORE GIVING UP.
C --------------------------------------------------------------------
 60   CONTINUE
      
      IF (NHCUT.EQ.10) THEN
C        HAVE REDUCED H TEN TIMES
         WRITE (LOUT,9100)
         GO TO 70
      ENDIF

      NHCUT = NHCUT + 1
      HMIN = 0.1D+0*HMIN
      H = 0.1D+0*H
      JSTART = -1  
      GO TO 20
C
 70   IF(DABS(T-TOUT).GT.1000.0D+0*UROUND) THEN
         DO 80 I = 1,N
            Y0(I) = Y(I,1)
 80      CONTINUE
         T0 = T

      ELSE
C        HAVE PASSED TOUT SO INTERPOLATE
         CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0)
         T0 = TOUT
         IDID = KFLAG
      ENDIF
      HO = H
      IF(KFLAG.NE.0) IDID = KFLAG
      RETURN
C -------------------------- END OF SUBROUTINE OVDRIV -----------------
 9000 FORMAT (' WARNING..  T + H = T ON NEXT STEP.')
 9010 FORMAT (/,/,' KFLAG = -2 FROM INTEGRATOR AT T = ',E16.8,'  H =',
     +       E16.8,/,
     +       '  THE REQUESTED ERROR IS SMALLER THAN CAN BE HANDLED',/,/)
 9020 FORMAT (/,/,' INTEGRATION HALTED BY DRIVER AT T = ',E16.8,/,
     +     '  ERROR TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION',
     +       /)
 9030 FORMAT (/,/,' KFLAG = -3 FROM INTEGRATOR AT T = ',E16.8,/,
     +       '  CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED',/)
 9040 FORMAT (/,/,' ILLEGAL INPUT.. RTOL .LE. 0.',/,/)
 9045 FORMAT (/,/,' ILLEGAL INPUT.. ATOL .LE. 0.',/,/)
 9050 FORMAT (/,/,' ILLEGAL INPUT.. N .LE. 0',/,/)
 9060 FORMAT (/,/,' ILLEGAL INPUT.. (T0-TOUT)*H .GE. 0.',/,/)
 9070 FORMAT (/,/,' ILLEGAL INPUT.. IDID =',I5,/,/)
 9080 FORMAT (/,/,' IDID = -1 ON INPUT WITH (T-TOUT)*H .GE. 0.',/,
     +       ' T =',E16.8,'   TOUT =',E16.8,'   H =',E16.8,/,
     +       ' INTERPOLATION WAS DONE AS ON NORMAL RETURN.',/,
     +       ' DESIRED PARAMETER CHANGES WERE NOT MADE.')
 9090 FORMAT (/,/,' ILLEGAL INPUT.. METHOD FLAG, MF, = ',I6,/,
     +       '         ALLOWED VALUES ARE 25 OR 26',/)
 9100 FORMAT (/,/,' PROBLEM APPEARS UNSOLVABLE WITH GIVEN INPUT',/,
     +       '         HMIN REDUCED BY A FACTOR OF 1.0E10',/,/)
 9110 FORMAT (/,/,' ILLEGAL VALUE FOR ITOL',/,/)
 9120 FORMAT (/,/,' ILLEGAL VALUE FOR MAXDER',/,/)
 9130 FORMAT (/,/,' NUMBER OF STEPS EXCEEDS MAXIMUM',/,/)
 9150 FORMAT (/,/,'IERR IS NON-ZERO BECAUSE OF AN ILLEGAL FUNCTION
     +CALL')
 9160 FORMAT (/,/,'STEPSIZE IS TOO SMALL')

      END 
C-----------------------------------------------------------------------------
      SUBROUTINE INTERP(N,JSTART,H,T,Y,TOUT,Y0)
  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C     .. SCALAR ARGUMENTS ..
      INTEGER JSTART,N
C     ..
C     .. ARRAY ARGUMENTS ..
      DIMENSION  Y(N,12),Y0(N)
C     ..
C     .. LOCAL SCALARS ..
      INTEGER I,J,L 
C     ..
C     .. INTRINSIC FUNCTIONS ..
C     ..
      DO 10 I = 1,N 
         Y0(I) = Y(I,1)
 10   CONTINUE
      L = JSTART + 2
      S = (TOUT-T)/H
      S1 = 1.0D+0
      DO 30 J = 2,L 
         S1 = S1* (S+DBLE(FLOAT(J-2)))/DBLE(FLOAT(J-1))
         DO 20 I = 1,N
            Y0(I) = Y0(I) + S1*Y(I,J)
 20      CONTINUE
 30   CONTINUE
      RETURN
C -------------- END OF SUBROUTINE INTERP --------------------------- 
      END 
C
      SUBROUTINE COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD)
      
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C --------------------------------------------------------------------
C     COSET IS CALLED BY THE INTEGRATOR AND SETS THE COEFFICIENTS USED
C     BY THE CONVENTIONAL BACKWARD DIFFERENTIATION SCHEME AND THE
C     MODIFIED EXTENDED BACKWARD DIFFERENTIATION SCHEME.  THE VECTOR
C     EL OF LENGTH NQ+1 DETERMINES THE BASIC BDF METHOD WHILE THE VECTOR
C     ELST OF LENGTH NQ+2 DETERMINES THE MEBDF.  THE VECTOR TQ OF
C     LENGTH 4 IS INVOLVED IN ADJUSTING THE STEPSIZE IN RELATION TO THE
C     TRUNCATION ERROR.  ITS VALUES ARE GIVEN BY THE PERTST ARRAY.  THE
C     VECTORS EL AND TQ BOTH DEPEND ON METH AND NQ.  THE  
C     COEFFICIENTS IN PERTST NEED TO BE GIVEN TO ONLY ABOUT ONE PERCENT
C     ACCURACY.  THE ORDER IN WHICH THE GROUPS APPEAR BELOW IS: 
C     COEFFICIENTS FOR ORDER NQ-1, COEFFICIENTS FOR ORDER NQ,
C     COEFFICIENTS FOR ORDER NQ+1. 
C ------------------------------------------------------------------- 
C     .. SCALAR ARGUMENTS ..
      INTEGER NQ
C     ..
C     .. ARRAY ARGUMENTS ..
      DIMENSION  EL(10),ELST(10),TQ(5)
C     ..
C     .. LOCAL SCALARS ..
      INTEGER K
C     ..
C     .. LOCAL ARRAYS ..
      DIMENSION  PERTST(8,3)
C     ..
C     .. INTRINSIC FUNCTIONS ..
C     ..
C     .. COMMON BLOCKS ..
      INTEGER MAXORD,NBSOL,NCOSET,NDEC,NFE,NJE,NPSET,NQUSED,NSTEP
C     ..
C     .. DATA STATEMENTS ..
      DATA PERTST(1,1)/1./,PERTST(2,1)/2./,PERTST(3,1)/4.5/,
     +     PERTST(4,1)/7.333/,PERTST(5,1)/10.42/,
     +     PERTST(6,1)/13.7/,PERTST(7,1)/17.15/,
     +     PERTST(8,1)/20.74/ 
      DATA PERTST(1,2)/2./,PERTST(2,2)/4.5/,
     +     PERTST(3,2)/7.333/,PERTST(4,2)/10.42/,
     +     PERTST(5,2)/13.7/,PERTST(6,2)/17.15/,
     +     PERTST(7,2)/20.74/,PERTST(8,2)/24.46/ 
      DATA PERTST(1,3)/4.5/,PERTST(2,3)/7.333/,
     +     PERTST(3,3)/10.42/,PERTST(4,3)/13.7/,
     +     PERTST(5,3)/17.15/,PERTST(6,3)/20.74/,
     +     PERTST(7,3)/24.46/,PERTST(8,3)/1./
C     ..
C ------------------------------------------------------------------- 
C     THE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO MACHINE ACCURACY.
C     THEIR DERIVATION IS GIVEN IN REFERENCE 1.
C ------------------------------------------------------------------- 
      IF (NQ.GT.MAXORD) MAXORD = NQ
      NCOSET = NCOSET + 1
      GO TO (10,20,30,40,50,60,70) NQ
  
 10   EL(1) = 1.0D+0
      ELST(1) = 1.5D+0
      ELST(3) = -0.5D+0
      GO TO 80
      
 20   EL(1) = 6.6666666666667D-01
      EL(3) = 3.3333333333333D-01
      ELST(1) = 9.5652173913043D-01
      ELST(3) = 2.1739130434782D-01
      ELST(4) = -1.7391304347826D-01
      GO TO 80
      
 30   EL(1) = 5.4545454545455D-01
      EL(3) = 4.5454545454545D-01
      EL(4) = 1.8181818181818D-01
      ELST(1) = 7.6142131979695D-01
      ELST(3) = 3.2994923857868D-01
      ELST(4) = 8.6294416243654D-02
      ELST(5) = -9.1370558375634D-02
      GO TO 80
      
 40   EL(1) = 0.48D+0
      EL(3) = 0.52D+0
      EL(4) = 0.28D+0
      EL(5) = 0.12D+0
      ELST(1) = 6.5733706517393D-01
      ELST(3) = 4.0023990403838D-01
      ELST(4) = 1.5793682526989D-01
      ELST(5) = 4.4382247101159D-02
      ELST(6) = -5.7576969212315D-02
      GO TO 80
      
 50   EL(1) = 4.3795620437956D-01
      EL(3) = 5.62043795620436D-01
      EL(4) = 3.43065693430656D-01
      EL(5) = 1.97080291970802D-01
      EL(6) = 8.75912408759123D-02
      ELST(1) = 5.9119243917152D-01
      ELST(3) = 4.4902473356122D-01
      ELST(4) = 2.1375427307460D-01
      ELST(5) = 9.0421610027481503D-02
      ELST(6) = 2.6409276761177D-02
      ELST(7) = -4.0217172732757D-02
      GO TO 80
  
 60   EL(1) = 4.08163265306120D-01
      EL(3) = 5.91836734693874D-01
      EL(4) = 3.87755102040813D-01
      EL(5) = 2.51700680272107D-01
      EL(6) = 1.49659863945577D-01
      EL(7) = 6.80272108843534D-02
      ELST(1) = 5.4475876041119D-01
      ELST(3) = 4.8525549636077D-01
      ELST(4) = 2.5789750131312D-01
      ELST(5) = 1.3133738525800D-01
      ELST(6) = 5.7677396763462D-02
      ELST(7) = 1.7258197643881D-02
      ELST(8) = -3.0014256771967D-02
      GO TO 80
  
 70   EL(1) = 3.85674931129476D-01
      EL(3) = 6.14325068870521D-01
      EL(4) = 4.21487603305783D-01
      EL(5) = 2.9292929292929D-01
      EL(6) = 1.96510560146923D-01
      EL(7) = 1.19375573921028D-01
      EL(8) = 5.50964187327820D-02
      ELST(1) = 5.0999746293734D-01
      ELST(3) = 5.1345839935281D-01
      ELST(4) = 2.9364346131937D-01
      ELST(5) = 1.6664672120553D-01
      ELST(6) = 8.8013735242353D-02
      ELST(7) = 3.9571794884069D-02
      ELST(8) = 1.2039080338722D-02
      ELST(9) = -2.3455862290154D-02
 80   DO 90 K = 1,3 
         TQ(K) = PERTST(NQ,K)
 90   CONTINUE
      TQ(4) = 0.5D+0*TQ(2)/DBLE(FLOAT(NQ))
      IF(NQ.NE.1) TQ(5)=PERTST(NQ-1,1)
      RETURN
C --------------------- END OF SUBROUTINE COSET --------------------- 
      END 
C
      SUBROUTINE PSET(Y,N,H,T,UROUND,EPSJAC,CON,MITER,IER,F,PDERV,
     +     NRENEW,YMAX,SAVE1,SAVE2,PW,PWCOPY,WRKSPC,ITOL,RTOL,ATOL,
     +     NPSET,NJE,NFE,NDEC,RSP,IA,JA,P,IP,NSP,ISP,NONZ,SAVF,NGRP,
     +     JGP,IGP,LOUT,IPAR,RPAR,IERR)
  
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C ------------------------------------------------------------------- 
C     PSET IS CALLED BY STIFF TO COMPUTE AND PROCESS THE MATRIX
C     I/(H*EL(1)) - J  WHERE J IS AN APPROXIMATION TO THE RELEVANT JACOBIAN
C     AND I IS THE IDENTITY MATRIX.  THIS MATRIX IS THEN SUBJECTED TO LU 
C     DECOMPOSITION IN PREPARATION FOR LATER SOLUTION OF LINEAR SYSTEMS
C     OF ALGEBRAIC EQUATIONS WITH LU AS THE COEFFICIENT MATRIX.  
C     IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION WITH
C     PSET USES THE FOLLOWING ..
C     EPSJAC = DSQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS.
C ******************************************************************* 
C     THE ARGUMENT NRENEW IS USED TO SIGNAL WHETHER OR NOT
C     WE REQUIRE A NEW JACOBIAN TO BE CALCULATED. 
C
C        IF NRENEW > 0 THEN WE REQUIRE A NEW J TO BE COMPUTED
C                  = 0 THEN USE A COPY OF THE LAST J COMPUTED
C ******************************************************************* 
C    THE FLAG IER IS SET TO 1 WHEN A STRUCTURAL SINGULARITY OCCURS 
C    AND TO 2 WHEN A NUMERICAL SINGULARITY OCCURS
C ******************************************************************* 
C 
C     .. SCALAR ARGUMENTS ..
      INTEGER IER,MITER,N,NRENEW
C     ..
C     .. ARRAY ARGUMENTS ..
      DIMENSION  PW(*),PWCOPY(*),SAVE2(N),WRKSPC(N),Y(N,12),
     &     YMAX(N),SAVE1(N),RTOL(1),ATOL(1),RSP(NSP),FTEM(1),SAVF(N),
     &     IPAR(*),RPAR(*)
      INTEGER P(N),IP(N),IA(N+1),JA(NONZ),ISP(NSP),IGP(NGRP+1),JGP(N)
C 
      LOGICAL ORDER
C     .. EXTERNAL SUBROUTINES ..
      EXTERNAL F,PDERV,ODRV,CDRV
C     ..
C     .. INTRINSIC FUNCTIONS ..
      INTRINSIC DABS,DMAX1,DSQRT
C     ..
C     .. COMMON BLOCKS ..
      COMMON /ORDERING/ORDER
      INTEGER MAXORD,NBSOL,NCOSET,NDEC,NFE,NJE,NPSET,NQUSED,NSTEP
c     
c   
      NPSET = NPSET + 1       
      IF (NRENEW.EQ.0) THEN
         IF (MITER.EQ.5 .OR. MITER .EQ. 6) THEN             
            KK=1
            KMIN = IA(1)  
            DO 275 J=1,N
               KMAX=IA(1+J)-1
               DO 270 K=KMIN,KMAX
                  I =JA(K)
                  PIJ=-PWCOPY(K)             
                  IF ( I .EQ. J) PIJ = PIJ +1.0D0/CON
                  PW(K) = PIJ
 270           CONTINUE
               KMIN=KMAX+1
 275        CONTINUE
         ENDIF
         GO TO 131
      ENDIF
      NJE = NJE + 1
 70   CONTINUE
C           
C.... FORM THE ITERATION MATRIX P =I/CON - J
C.... MITER = 5, CALL PDERV ROUTINE AND FORM P
C...
      IF (MITER .EQ. 5) THEN           
         KMIN = IA(1)
         DO 130 J = 1,N
            KMAX = IA(1+J) -1
            DO 110 I = 1,N
 110           PWCOPY(I) = 0.0D0
            CALL PDERV(N,T,Y,J,PWCOPY,IPAR,RPAR,IERR)
            DO 120 K = KMIN,KMAX
               II = JA(K)
               PW(K) = -PWCOPY(II)                   
               IF (II .EQ. J) PW(K) = PW(K) + 1.0D0/CON 
 120        CONTINUE
            KMIN = KMAX + 1           
 130     CONTINUE
      ELSE 
C....
C.... NUMERICAL APRROXIMATION TO THE SPARSE JACOBIAN, MAKE NGRP CALS 
C.... TO F TO APPROXIMATE THE JACOBIAN AND FORM THE MATRIX P
C....         
         DO 35 I=1,N
            IF(ITOL.EQ.2) THEN
               YMAX(I)=dabs(Y(I,1))*RTOL(1)+ATOL(1)
            ELSE IF (ITOL.EQ.3) THEN
               YMAX(I)=dabs(Y(I,1))*RTOL(1)+ATOL(I)
            ELSE IF(ITOL.EQ.4) THEN
               YMAX(I)=dabs(Y(I,1))*RTOL(I)+ATOL(1)
            ELSE IF(ITOL.EQ.5) THEN
               YMAX(I)=dabs(Y(I,1))*RTOL(I)+ATOL(I)
            END IF
 35      CONTINUE
         D = 0.D0 
         DO 41 I = 1,N
            D = D + (SAVE2(I)*YMAX(I))**2
 41      CONTINUE
         IF (ITOL .EQ. 1) D= D*RTOL(1)**2
         D = DSQRT(D/DFLOAT(N))
         R0 = DABS(H)*D*(1.0D+03)*UROUND*DFLOAT(N)
         IF (R0.EQ.0.D0) R0 = 1.D0 
         JMIN = IGP(1)           
         DO 240 NG =1,NGRP
            JMAX= IGP(1+NG)-1              
            DO 210 J=JMIN, JMAX
               JJ= JGP(J)                     
               SAVE1(JJ) =Y(JJ,1)
               YI = Y(JJ,1)                 
               R = DMAX1(EPSJAC*DABS(YI),R0/YMAX(JJ)) 
 210           Y(JJ,1)=Y(JJ,1)+R                  
            CALL F(N,T,Y,WRKSPC,IPAR,RPAR,IERR)                  
            DO 230 J=JMIN,JMAX
               JJ = JGP(J)               
               Y(JJ,1)=SAVE1(JJ)
               YJJ = Y(JJ,1)
               R = DMAX1(EPSJAC*DABS(YJJ),R0/YMAX(JJ))
               D = CON/R
               K1 = IA(JJ)
               K2 = IA(JJ+1) -1                     
               DO 220 K=K1,K2
                  I = JA(K)
                  TEMPRY = WRKSPC(I) - SAVE2(I)
                  PWCOPY(k)=TEMPRY/R
                  PW(k) = -PWCOPY(k)
                  IF (I .EQ. JJ ) PW(K) = PW(K) + 1.0D0/CON
 220           CONTINUE
 230        CONTINUE
            JMIN = JMAX + 1
 240     CONTINUE      
         NFE = NFE + NGRP   
      ENDIF        
 131  CONTINUE
C....      
C.... NUMERICAL FACTORIZATION OF THE ITERATION MATRIX P
C....        
      DO 132 I=1,N
 132     WRKSPC(I) = 0.0D0
      CALL CDRV(N,P,P,IP,IA,JA,PW,WRKSPC,WRKSPC,NSP,ISP,
     +        RSP,IEXTRA,2,IFLAG)                 
C     
      IF (IFLAG .NE. 0) THEN           
         WRITE(LOUT,1333) IFLAG
 1333    FORMAT('FATAL ERROR IN THE CDRV ROUTINE (SPARSE
     +        NUMERICAL FACTORIZATION) IFLAG RETURN IS ',I3)
         IER = 1
         ORDER = .FALSE.
      ENDIF                                     
      NDEC = NDEC + 1
C     
      RETURN
      END 
C ---------------------- END OF SUBROUTINE PSET --------------------- 

      SUBROUTINE ERRORS(N,TQ,EDN,E,EUP,BND,EDDN)  

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C     ***************************************************
C
C     THIS ROUTINE CALCULATES ERRORS USED IN TESTS
C     IN STIFF .
C
C     ***************************************************
C     .. SCALAR ARGUMENTS ..
      INTEGER N
C     ..
C     .. ARRAY ARGUMENTS ..
      DIMENSION  TQ(5)
C     ..
C     .. LOCAL SCALARS ..
C     ..
C     .. INTRINSIC FUNCTIONS ..
C     ..
      SQHOL = DBLE(FLOAT(N))
      EDN = TQ(1)*TQ(1)*SQHOL 
C
C     ** ERROR ASSOCIATED WITH  METHOD OF ORDER ONE LOWER.
C
      E = TQ(2)*TQ(2)*SQHOL
C
C     ** ERROR ASSOCIATED WITH PRESENT ORDER
C
      EUP = TQ(3)*TQ(3)*SQHOL 
C
C     ** ERROR ASSOCIATED WITH HIGHER ORDER METHOD
C
      BND = TQ(4)*TQ(4)*SQHOL*0.5D+0

      EDDN=TQ(5)*TQ(5)*SQHOL
C
C     ** ERROR ASSOCIATED WITH METHOD OF ORDER TWO LOWER.
      RETURN
      END 
C-----------------------------------------------------------------------------

      SUBROUTINE PRDICT(T,H,Y,L,N,YPRIME,NFE,F,IPAR,RPAR,IERR)  

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C **********************************************************************
C     PREDICTS A VALUE FOR Y AT (T+H) GIVEN THE HISTORY ARRAY AT T 
C     THEN EVALUATES THE DERIVATIVE AT THIS POINT, THE RESULT OF THIS 
C     EVALUATION BEING STORED IN YPRIME()
C **********************************************************************
  
C     .. SCALAR ARGUMENTS ..
      INTEGER L,N
C     ..
C     .. ARRAY ARGUMENTS ..
      DIMENSION  Y(N,12),YPRIME(N),IPAR(*),RPAR(*)
C     ..
C     .. LOCAL SCALARS ..
      INTEGER I,J2
C     ..
C     .. EXTERNAL SUBROUTINES ..
      EXTERNAL F
C     ..
C     .. COMMON BLOCKS ..
      INTEGER NFE
C     ..      
      DO 20 J2 = 2,L
         DO 10 I = 1,N
            Y(I,1) = Y(I,1) + Y(I,J2)         
 10      CONTINUE
 20   CONTINUE
      T = T + H
      CALL F(N,T,Y,YPRIME,IPAR,RPAR,IERR)
      NFE = NFE + 1 
      RETURN
      END 

C-----------------------------------------------------------------------------

      SUBROUTINE ITRAT2(QQQ,Y,N,T,HBETA,ERRBND,ARH,CRATE,TCRATE,
     +     M,WORKED ,YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,MF,LMB,ITOL,
     +     RTOL,ATOL,HUSED,NBSOL,NFE,NQUSED,F,RSP,IA,JA,P,IP,NSP,
     +     ISP,NONZ,LOUT,IPAR,RPAR,IERR) 
C     
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C     .. SCALAR ARGUMENTS ..
      INTEGER M,N,ITOL
      LOGICAL WORKED
C     ..
C     .. ARRAY ARGUMENTS ..
      DIMENSION ARH(N),ERROR(N),PW(nonz),SAVE1(N),SAVE2(N),Y(N,12),
     &     YMAX(N),RTOL(1),ATOL(1),SCALE(N),RSP(NSP),RPAR(*)
      INTEGER P(N),IP(N),ISP(NSP),IA(N+1),JA(NONZ),IPAR(*)
      INTEGER MAXORD,NBSOL,NCOSET,NDEC,NFE,NJE,NPSET,NQUSED,NSTEP
C      
C     .. EXTERNAL SUBROUTINES ..
      EXTERNAL F,CDRV
C     ..
C     .. INTRINSIC FUNCTIONS ..
      INTRINSIC DMAX1,DMIN1      
C     ..
C     .. DATA STATEMENTS ..
      DATA  ZERO/0.0D+0/
C     ..
C            
      DO 5 I=1,N
         AYI = DABS(Y(I,1))
         IF(ITOL.EQ.1) THEN
            SCALE(I) = YMAX(I)
         ELSE IF(ITOL.EQ.2) THEN
            SCALE(I) = RTOL(1)*AYI + ATOL(1)
         ELSE IF(ITOL.EQ.3) THEN
            SCALE(I) = RTOL(1)*AYI + ATOL(I)
         ELSE IF(ITOL.EQ.4) THEN
            SCALE(I) = RTOL(I)*AYI + ATOL(1)
         ELSE IF(ITOL.EQ.5) THEN
            SCALE(I) = RTOL(I)*AYI + ATOL(I)
         ENDIF
 5    CONTINUE
      IF(LMB.EQ.1) GOTO 25
      
      DO 10 I = 1,N
         SAVE1(I) = (-SAVE1(I)+HBETA*SAVE2(I)+ARH(I))/QQQ
 10   CONTINUE
      IFLAG = 0      
      IF ( MF .EQ. 25 .OR. MF .EQ. 26 ) THEN
         CALL CDRV(N,P,P,IP,IA,JA,PW,SAVE1,SAVE1,NSP,ISP,
     +        RSP,IEXTRA,4,IFLAG)
         IF (IFLAG .NE. 0) THEN
            WRITE(LOUT,2222) IFLAG
 2222       FORMAT('FALTAL ERROR IN THE CDRV ROUTINE 
     +           (SPARSE SOLVER) IFLAG RETURN IS ' , I3)
            STOP
         ENDIF 
         NBSOL = NBSOL + 1
      ENDIF
      
      D = ZERO
      DO 20 I = 1,N 
         ERROR(I) = ERROR(I) + SAVE1(I)
         D = D + (SAVE1(I)/(SCALE(I)))**2
         SAVE1(I) = Y(I,1) + ERROR(I)
 20   CONTINUE
      IF(ITOL.EQ.1) D=D/(RTOL(1)**2)
      TCRATE = TCRATE + CRATE 
      D1 = D
      M = 1
      CALL F(N,T,SAVE1,SAVE2,IPAR,RPAR,IERR) 
      NFE = NFE + 1 

 25   CONTINUE
      WORKED = .TRUE.
 30   CONTINUE      
      DO 40 I = 1,N 
         SAVE1(I)=(-SAVE1(I)+HBETA*SAVE2(I)+ARH(I))/qqq        
 40   CONTINUE
C     
C     IF WE ARE HERE THEN PARTIALS ARE O.K.
C
      IFLAG  = 0
      IF( MF.EQ. 25 .OR. MF .EQ. 26) THEN          
         CALL CDRV(N,P,P,IP,IA,JA,PW,SAVE1,SAVE1,NSP,ISP,
     +        RSP,IEXTRA,4,IFLAG)
         IF (IFLAG .NE. 0) THEN
            WRITE(LOUT,3333) IFLAG
 3333       FORMAT('FALTAL ERROR IN THE CDRV ROUTINE 
     +           (SPARSE SOLVER) IFLAG RETURN IS ', I3)            
            STOP
         ENDIF          
         NBSOL=NBSOL + 1
      ENDIF
C
C     WE NOW CALCULATE A WEIGHTED RMS TYPE NORM
C
      D = ZERO
      DO 50 I = 1,N 
         ERROR(I) = ERROR(I) + SAVE1(I)
         D = D + (SAVE1(I)/(SCALE(I)))**2
         SAVE1(I) = Y(I,1) + ERROR(I)
 50   CONTINUE

      IF(ITOL.EQ.1) D=D/(RTOL(1)**2)
C ------------------------------------------------------------------- 
C     TEST FOR CONVERGENCE.  IF M.GT.0 , AN ESTIMATE OF THE CONVERGENCE
C     RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST. 
C -------------------------------------------------------------------       
      IF (M.NE.0) THEN
         IF (D1.NE.ZERO) CRATE = DMAX1(0.9D+0*CRATE,D/D1)
      ENDIF
      
      TCRATE = TCRATE + CRATE 
      IF ((D*DMIN1(1.0D+0,2.0D+0*CRATE)).LT.ERRBND/DFLOAT(NQUSED))
     +     RETURN
      IF (M.NE.0) THEN
         IF (D.GT.D1) THEN
            WORKED = .FALSE.
            RETURN
         ENDIF
      ENDIF
      
      D1 = D
      IF (M.EQ.4) THEN
         WORKED = .FALSE.
         RETURN
      ENDIF
      
      M = M + 1
      CALL F(N,T,SAVE1,SAVE2,IPAR,RPAR,IERR)
      NFE = NFE + 1
      GO TO 30
      END 
      
C-----------------------------------------------------------------------------

      SUBROUTINE STIFF(H,HMAX,HMIN,JSTART,KFLAG,MF,T,TOUT,TEND,Y,N,
     +     YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,PWCOPY,YHOLD,YNHOLD,ARH,
     +     LOUT,MAXDER,ITOL,RTOL,ATOL,F,PDERV,NQUSED,NSTEP,NFAIL,NFE,
     +     NJE,NDEC,NBSOL,NPSET,NCOSET,MAXORD,UROUND,EPSJAC,HUSED,RSP,
     +     IA,JA,P,IP,NSP,ISP,NONZ,SAVF,NGP,JGP,IGP,IPAR,RPAR,IERR)

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C     ------------------------------------------------------------------
C     THE SUBROUTINE STIFF PERFORMS ONE STEP OF THE INTEGRATION OF AN 
C     INITIAL VALUE PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL
C     EQUATIONS OR LINEARLY IMPLICIT DIFFERENTIAL ALGEBRAIC EQUATIONS.
C     COMMUNICATION WITH STIFF IS DONE WITH THE FOLLOWING VARIABLES.. 
C     Y      AN N BY LMAX+3 ARRAY CONTAINING THE DEPENDENT VARIABLES
C              AND THEIR BACKWARD DIFFERENCES.  MAXDER (=LMAX-1) IS THE 
C              MAXIMUM ORDER AVAILABLE.  SEE SUBROUTINE COSET.
C              Y(I,J+1) CONTAINS THE JTH BACKWARD DIFFERENCE OF Y(I)
C     T      THE INDEPENDENT VARIABLE. T IS UPDATED ON EACH STEP TAKEN.
C     H      THE STEPSIZE TO BE ATTEMPTED ON THE NEXT STEP. 
C              H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING
C              THE PROBLEM.  H CAN BE EITHER POSITIVE OR NEGATIVE BUT 
C              ITS SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM.
C     HMIN   THE MINIMUM AND MAXIMUM ABSOLUTE VALUE OF THE STEPSIZE
C     HMAX   TO BE USED FOR THE STEP.  THESE MAY BE CHANGED AT ANY
C              TIME BUT WILL NOT TAKE EFFECT UNTIL THE NEXT H CHANGE. 
C     RTOL,ATOL  THE ERROR BOUNDS. SEE DESCRIPTION IN OVDRIV.
C     UROUND THE UNIT ROUNDOFF OF THE MACHINE.
C     N      THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS.
C     MF     THE METHOD FLAG.  MUST BE SET TO 21,22,23 OR 24 AT PRESENT
C     KFLAG  A COMPLETION FLAG WITH THE FOLLOWING MEANINGS..
C                  0  THE STEP WAS SUCCESSFUL
C                 -1  THE REQUESTED ERROR COULD NOT BE ACHIEVED
C                       WITH ABS(H) = HMIN.
C                 -2  THE REQUESTED ERROR IS SMALLER THAN CAN
C                       BE HANDLED FOR THIS PROBLEM.
C                 -3  CORRECTOR CONVERGENCE COULD NOT BE
C                       ACHIEVED FOR ABS(H)=HMIN. 
C            ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF T AND 
C            THE Y ARRAY ARE AS AT THE BEGINNING OF THE LAST
C            STEP ATTEMPTED, AND H IS THE LAST STEP SIZE ATTEMPTED.
C     JSTART  AN INTEGER USED ON INPUT AND OUTPUT.
C          ON INPUT IT HAS THE FOLLOWING VALUES AND MEANINGS.. 
C              0  PERFORM THE FIRST STEP.
C          .GT.0  TAKE A NEW STEP CONTINUING FROM THE LAST
C          .LT.0  TAKE THE NEXT STEP WITH A NEW VALUE OF H OR N.
C          ON EXIT, JSTART IS NQUSED, THE ORDER OF THE METHOD LAST USED. 
C     YMAX     AN ARRAY OF N ELEMENTS WITH WHICH THE ESTIMATED LOCAL
C              ERRORS IN Y ARE COMPARED 
C     ERROR    AN ARRAY OF N ELEMENTS.
C     SAVE1,2  TWO ARRAYS OF WORKING SPACE BOTH OF LENGTH N.
C     PW       A BLOCK OF LOCATIONS USED FOR PARTIAL DERIVATIVES
C     IPIV     AN INTEGER ARRAY OF LENGTH N USED FOR PIVOT INFORMATION.
C     JNEWIM   IS TO INDICATE IF PRESENT ITERATION MATRIX
C                WAS FORMED USING A NEW J OR OLD J.
C     JSNOLD   KEEPS TRACK OF NO. OF STEPS TAKEN WITH
C                PRESENT ITERATION MATRIX (BE IT FORMED BY
C                A NEW J OR NOT).
C     AVNEWJ   STORES VALUE FOR AVERAGE CRATE WHEN ITERATION
C                MATRIX WAS FORMED BY A NEW J.
C     AVOLDJ   STORES VALUE FOR AVERAGE CRATE WHEN ITERATION
C                MATRIX WAS FORMED BY AN OLD J.
C     NRENEW   FLAG THAT IS USED IN COMMUNICATION WITH SUBROUTINE PSET.
C                IF  NRENEW > 0  THEN FORM A NEW JACOBIAN BEFORE
C                                COMPUTING THE COEFFICIENT MATRIX FOR 
C                                THE NEWTON-RAPHSON ITERATION
C                           = 0  FORM THE COEFFICIENT MATRIX USING A
C                                COPY OF AN OLD JACOBIAN
C     NEWPAR   FLAG USED IN THIS SUBROUTINE TO INDICATE IF A JACOBIAN 
C              HAS BEEN EVALUATED FOR THE CURRENT STEP
C **********************************************************************  
C     .. SCALAR ARGUMENTS ..
      INTEGER JSTART,KFLAG,LOUT,MF,N,ITOL,NZ
C     ..
C     .. ARRAY ARGUMENTS ..
      DIMENSION HSTPSZ(2,14)
      DIMENSION  ARH(N),ERROR(N),PW(*),PWCOPY(*),SAVE1(N),SAVE2(N),
     +     Y(N,12), YHOLD(N,12),YMAX(N),YNHOLD(N,2),RTOL(1),ATOL(1),
     +     SCALE(N),RSP(NSP),IPAR(*),RPAR(*),SAVF(N) 
      INTEGER P(N),IP(N),IA(N+1),JA(NONZ),ISP(NSP),JGP(N),IGP(NONZ)
C     ..
C     .. LOCAL SCALARS ..
      INTEGER I,IER,IITER,IITER2,IREDO,ITST,J,J1,J2,KFAIL,LL,LMAX,
     +     M3STEP,MAXDER,METH,MFOLD,MITER,NEWQ, NSP
C     
      LOGICAL FINISH,OVRIDE,WORKED
C     ..
C     .. LOCAL ARRAYS ..
      DIMENSION  EL(10),ELST(10),TQ(5)
C     ..
C     .. EXTERNAL SUBROUTINES ..
      EXTERNAL COSET,CPYARY,ERRORS,F,PDERV,HCHOSE,ITRAT2,
     +     PRDICT,PSET,RSCALE,CDRV
C     ..
C     .. INTRINSIC FUNCTIONS ..
      INTRINSIC DABS,DMAX1,DMIN1
C     ..
C     .. COMMON BLOCKS ..
      COMMON/STPSZE/HSTPSZ
      INTEGER IDOUB,ISAMP,IWEVAL,JCHANG,JSINUP,JSNOLD,L,M1,M2,
     +     MAXORD, MEQC1,MEQC2,MQ1TMP,MQ2TMP,NBSOL,NCOSET,
     +     NDEC,NEWPAR,NFE,NJE,NPSET,NQ,NQUSED,NRENEW,NSTEP,NT

      LOGICAL CFAIL,JNEWIM,SAMPLE
C     ..
C     .. SAVE STATEMENT ..
C     ..
      SAVE  MFOLD,HOLD,LMAX,EDN,EUP,BND,EDDN,EL,TQ,ELST,E,IWEVAL,
     +     NEWPAR,NRENEW,JSINUP,JSNOLD,IDOUB,JCHANG,L,NQ,MEQC1,QI,QQQ,
     +     MEQC2,MQ1TMP,MQ2TMP,ISAMP,RH,RMAX,TOLD,CRATE1,CRATE2,IER,
     +     TCRAT1,TCRAT2,AVNEW2,AVOLD2,AVNEWJ,AVOLDJ,MITER,IBND,UPBND,
     +     CFAIL,KFAIL,RC,JNEWIM,VTOL,SAMPLE,IEMB,OLDLO,IREDO,OVRIDE


C     .. DATA STATEMENTS ..
      
      DATA  EL(2),ELST(2),OLDLO/3*1.0D+0/
      DATA  ZERO,ONE/0.0D+0,1.0D+0/
C     ..
      TOLD = T
      KFLAG = 0
      IF (JSTART.GT.0) GO TO 60
      IF (JSTART.NE.0) GO TO 30
C     ------------------------------------------------------------------
C     ON THE FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL YDOT
C     IS CALCULATED.  RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE
C     INCREASED IN A SINGLE STEP.  RMAX IS SET EQUAL TO 1.D4 INITIALLY
C     TO COMPENSATE FOR THE SMALL INITIAL H, BUT THEN IS NORMALLY = 10.
C     IF A FAILURE OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST),
C     RMAX IS SET AT 2. FOR THE NEXT INCREASE.
C     ------------------------------------------------------------------
      CALL F(N,T,Y,SAVE1,IPAR,RPAR,IERR)
      IF (IERR .NE. 0) GOTO 8000
      DO 10 I = 1,N 
        Y(I,2) = H*SAVE1(I)
 10   CONTINUE
      METH = 2
      MITER = MF - 10*METH
      IBND=5
      UPBND=0.1D+0
      NQ = 1
      NQUSED=NQ
      IER=0
      L = 2
      IDOUB = 3
      KFAIL = 0
      RMAX = 10000.0D+0
      RC = ZERO
      CRATE1 = 0.1D+0
      CRATE2 = 0.1D+0
      JSNOLD = 0
      JNEWIM = .TRUE.
      TCRAT1 = ZERO 
      TCRAT2 = ZERO 
      VTOL=DMAX1(RTOL(1),ATOL(1))/10.0D+0
      DO 15 I=1,12
         HSTPSZ(1,I)=1.0D+0 
         HSTPSZ(2,I)=VTOL
 15   CONTINUE
      HOLD = H
      MFOLD = MF
      NSTEP = 0
      NFE = 1
      NJE = 0
      NDEC = 0
      NPSET = 0
      NCOSET = 0
      MAXORD = 1
      NFAIL = 0
      CFAIL = .TRUE.
      AVNEWJ = ZERO 
      AVOLDJ = ZERO 
      AVNEW2 = ZERO 
      AVOLD2 = ZERO 
      SAMPLE = .FALSE.
      ISAMP = 0
      IEMB=0
C     **************************************************
C     CFAIL=.TRUE. ENSURES THAT WE CALCULATE A NEW
C     J ON THE FIRST CALL.
C     **************************************************
      MEQC1 = 0
      MEQC2 = 0
      MQ1TMP = 0
      MQ2TMP = 0
      NBSOL = 0
      HUSED = H
C     -----------------------------------------------------------------
C     IF THE CALLER HAS CHANGED N , THE CONSTANTS E, EDN, EUP
C     AND BND MUST BE RESET.  E IS A COMPARISON FOR ERRORS AT THE
C     CURRENT ORDER NQ.  EUP IS TO TEST FOR INCREASING THE ORDER,
C     EDN FOR DECREASING THE ORDER.  BND IS USED TO TEST FOR CONVERGENCE
C     OF THE CORRECTOR ITERATES.   IF THE CALLER HAS CHANGED H, Y MUST 
C     BE RE-SCALED.  IF H IS CHANGED, IDOUB IS SET TO L+1 TO PREVENT
C     FURTHER CHANGES IN H FOR THAT MANY STEPS.
C     -----------------------------------------------------------------
      CALL COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD)
      LMAX = MAXDER + 1
      RC = RC*EL(1)/OLDLO
      OLDLO = EL(1) 
      IWEVAL = MITER
      NRENEW = 1
      NEWPAR = 0
C     ***************************************************** 
C     NRENEW AND NEWPAR ARE TO INSTRUCT ROUTINE THAT
C     WE WISH A NEW J TO BE CALCULATED FOR THIS STEP.
C     ***************************************************** 
      CALL ERRORS(N,TQ,EDN,E,EUP,BND,EDDN)
      DO 20 I = 1,N 
         ARH(I) = EL(2)*Y(I,1)
 20   CONTINUE
      CALL CPYARY(N*L,Y,YHOLD)
      QI = H*EL(1)
      QQ = ONE/QI
      CALL PRDICT(T,H,Y,L,N,SAVE2,NFE,F,IPAR,RPAR,IERR)      
      IF (IERR .NE. 0 ) GOTO 8000
      GO TO 110
C     >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C     DIFFERENT PARAMETERS ON THIS CALL        <
C     <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 30   CALL CPYARY(N*L,YHOLD,Y)
      IF (MF.NE.MFOLD) THEN
         METH = MF/10
         MITER = MF - 10*METH 
         MFOLD = MF 
         IWEVAL = MITER
      ENDIF
      IF(NSTEP.GT.0) GOTO 35
      NJE = 0
      NFE = 1
      CFAIL = .TRUE.
      NEWPAR = 0
      MQ1TMP = 0
      MQ2TMP = 0
      MEQC1 = 0
      MEQC2 = 0
      TCRAT1 = 0.0D+0
      TCRAT2 = 0.0D+0
      CRATE1 = 1.0D-1
      CRATE2 = 1.0D-1
      NSTEP = 0
      NBSOL = 0
      NPSET = 0
      NCOSET = 0
      NDEC = 0
 35   CONTINUE

      IF (H.NE.HOLD) THEN
         RH = H/HOLD
         H = HOLD
         IREDO = 3
         GO TO 50
      ELSE
         GO TO 60
      ENDIF
C     *********************************************
C     RE-SCALE Y AFTER A CHANGE OF STEPSIZE   *
C     *********************************************
 40   RH = DMAX1(RH,HMIN/DABS(H))
 50   RH = DMIN1(RH,HMAX/DABS(H),RMAX)
      CALL RSCALE(N,L,RH,Y)
      RMAX = 10.0D+0
      JCHANG = 1
      H = H*RH
      RC = RC*RH
      IF (JSNOLD.GT.IBND) THEN
         CFAIL = .TRUE.
         NEWPAR = 0 
         RC = ZERO
C **********************************************************************
C        CFAIL=TRUE AND NEWPAR=0 SHOULD FORCE A NEW J TO BE EVALUATED 
C        AFTER 7 STEPS WITH AN OLD J, IF WE HAVE HAD A FAILURE OF ANY
C        KIND ON THE FIRST, SECOND OR THIRD STAGE OF THE CURRENT STEP 
C **********************************************************************
      ENDIF
      IDOUB = L + 1 
      CALL CPYARY(N*L,Y,YHOLD)
  
 60   IF (DABS(RC-ONE).GT.UPBND) IWEVAL = MITER
      HUSED = H
C     ------------------------------------------------------------------
C     THIS SECTION COMPUTES THE PREDICTED VALUES OF Y
C     AND THE RHS, ARH, FOR USE IN THE NEWTON ITERATION SCHEME.
C     RC IS THE RATIO OF THE NEW TO OLD VALUES OF THE COEFFICIENT
C     H*EL(1). WHEN RC DIFFERS FROM 1 BY MORE THAN 20 PERCENT, IWEVAL IS
C     SET TO MITER TO FORCE THE PARTIALS TO BE UPDATED.
C     ------------------------------------------------------------------
      QI = H*EL(1)
      QQ = ONE/QI
      DO 70 I = 1,N 
         ARH(I) = EL(2)*Y(I,1) 
 70   CONTINUE
      DO 90 J1 = 2,NQ
         JP1 =J1+1
         DO 80 I = 1,N
            ARH(I) = ARH(I) + EL(JP1)*Y(I,J1)
 80      CONTINUE
 90   CONTINUE
      
      IF (JCHANG.EQ.1) THEN
C        IF WE HAVE CHANGED STEPSIZE THEN PREDICT A VALUE FOR Y(T+H)
C        AND EVALUATE THE DERIVATIVE THERE (STORED IN SAVE2())
         CALL PRDICT(T,H,Y,L,N,SAVE2,NFE,F,IPAR,RPAR,IERR)         
         IF (IERR .NE. 0 ) GOTO 8000
      ELSE
C        ELSE USE THE VALUES COMPUTED FOR THE SECOND BDF FROM THE LAST
C        STEP. Y( ,LMAX+2) HOLDS THE VALUE FOR THE DERIVATIVE AT (T+H)
C        AND Y( ,LMAX+3) HOLDS THE APPROXIMATION TO Y AT THIS POINT.         
         LMP2=LMAX+2
         LMP3=LMAX+3
         DO 100 I = 1,N
            SAVE2(I) = Y(I,LMP2)
            Y(I,1) = Y(I,LMP3)
 100     CONTINUE
         T = T + H
      ENDIF
      
 110  IF (IWEVAL.LE.0) GO TO 120
C ------------------------------------------------------------------- 
C     IF INDICATED, THE MATRIX P = I/(H*EL(2)) - J IS RE-EVALUATED BEFORE
C     STARTING THE CORRECTOR ITERATION.  IWEVAL IS SET = 0 TO INDICATE
C     THAT THIS HAS BEEN DONE. P IS COMPUTED AND PROCESSED IN PSET.
C     THE PROCESSED MATRIX IS STORED IN PW
C ------------------------------------------------------------------- 
      IWEVAL = 0
      RC = ONE
      IITER = MEQC1 - MQ1TMP
      IITER2 = MEQC2 - MQ2TMP 
      IF (JNEWIM) THEN
         IF (JSNOLD.GE.3) THEN
            AVNEWJ = TCRAT1/DBLE(FLOAT(IITER))
            AVNEW2 = TCRAT2/DBLE(FLOAT(IITER2))
         ELSE
            AVNEWJ = ONE
            AVNEW2 = ONE
         ENDIF
      ELSE
C
C          MATRIX P WAS FORMED WITH A COPY OF J
C
         IF (JSNOLD.GE.3) THEN
            AVOLDJ = TCRAT1/DBLE(FLOAT(IITER))
            AVOLD2 = TCRAT2/DBLE(FLOAT(IITER2))
            IF (AVOLDJ.LT.AVNEWJ) THEN
               AVNEWJ = AVOLDJ
               
            ELSE IF (((DABS(AVOLDJ-AVNEWJ)).GT.0.3D+0) .OR. 
     +          ((AVOLDJ.GT.0.85D+0).AND. (AVOLDJ.NE.ONE))) THEN
C
C              SINCE IN CERTAIN INSTANCES AVOLDJ WILL
C              BE 1.0 AND THERE WILL BE NO NEED TO
C              UPDATE J.
C
               CFAIL = .TRUE. 
               CRATE1 = 0.1D+0
               CRATE2 = 0.1D+0
            ENDIF
         ELSE
            CFAIL = .TRUE.
            CRATE1 = 0.1D+0
            CRATE2 = 0.1D+0
C
C           ********************************************************* 
C           IF WE HAVE REACHED HERE THINGS MUST HAVE GONE WRONG
C           ********************************************************* 
C
         ENDIF
      ENDIF
  
      TCRAT1 = ZERO 
      TCRAT2 = ZERO 
      IF (CFAIL) THEN
         NRENEW = 1
         NEWPAR = 1
         JSINUP = -1
         JNEWIM = .TRUE.
      ENDIF
      
      CFAIL = .FALSE.
      JSNOLD = 0
      MQ1TMP = MEQC1
      MQ2TMP = MEQC2
      
      CALL PSET(Y,N,H,T,UROUND,EPSJAC,QI,MITER,IER,F,PDERV,NRENEW,YMAX,
     +     SAVE1,SAVE2,PW,PWCOPY,ERROR,ITOL,RTOL,ATOL,NPSET,NJE,NFE,
     +     NDEC,RSP,IA,JA,P,IP,NSP,ISP,NONZ,SAVF,NGP,JGP,IGP,LOUT,
     +     IPAR,RPAR,IERR)
C     
      IF (IERR .NE. 0 ) GOTO 8000
      QQQ=QI
C     NOTE THAT ERROR() IS JUST BEING USED AS A WORKSPACE BY PSET
      IF (IER.NE.0) THEN
C     IF IER>0 THEN WE HAVE HAD A SINGULARITY IN THE ITERATION MATRIX
         IJUS=1
         RED=0.5D+0
         NFAIL = NFAIL + 1
         GO TO 450
      ENDIF
      
 120  DO 130 I = 1,N
         SAVE1(I) = Y(I,1)
         ERROR(I) = ZERO
 130  CONTINUE
      M1 = 0
C **********************************************************************
C     UP TO 4 CORRECTOR ITERATIONS ARE TAKEN.  A CONVERGENCE TEST IS MADE
C     ON THE R.M.S. NORM OF EACH CORRECTION ,USING BND, WHICH DEPENDS
C     ON ATOL AND RTOL.  THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE
C     VECTOR  ERROR(I).  THE Y ARRAY IS NOT ALTERED IN THE CORRECTOR
C     LOOP. THE UPDATED Y VECTOR IS STORED TEMPORARILY IN SAVE1.
C **********************************************************************
      IF (.NOT.SAMPLE) THEN
         CALL ITRAT2(QQQ,Y,N,T,QI,BND,ARH,CRATE1,TCRAT1,M1,
     +        WORKED,YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,MF,
     +        1,ITOL,RTOL,ATOL,HUSED,NBSOL,NFE,NQUSED,F,
     +        RSP,IA,JA,P,IP,NSP,ISP,NONZ,LOUT,IPAR,RPAR,IERR)
         IF (IERR .NE. 0 ) GOTO 8000
         ITST = 2 
      ELSE
         CALL ITRAT2(QQQ,Y,N,T,QI,BND,ARH,CRATE1,TCRAT1,M1,
     +        WORKED,YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,MF,
     +        0,ITOL,RTOL,ATOL,HUSED,NBSOL,NFE,NQUSED,F,
     +        RSP,IA,JA,P,IP,NSP,ISP,NONZ,LOUT,IPAR,RPAR,IERR)
         IF (IERR .NE. 0 ) GOTO 8000
         ITST = 3         
      ENDIF
  
      MEQC1 = MEQC1 + M1 + 1
C
C       NOW TEST TO SEE IF IT WAS SUCCESSFUL OR NOT
C
      IF (.NOT.WORKED) THEN         
         NFAIL = NFAIL + 1
C **********************************************************************
C        THE CORRECTOR ITERATION FAILED TO CONVERGE IN 4 TRIES. IF
C        PARTIALS ARE NOT UP TO DATE, THEY ARE RE-EVALUATED FOR THE
C        NEXT TRY. OTHERWISE THE Y ARRAY IS REPLACED BY ITS VALUES
C        BEFORE PREDICTION AND H IS REDUCED IF POSSIBLE. IF NOT A
C        NON-CONVERGENCE EXIT IS TAKEN
C **********************************************************************
         IF (IWEVAL.EQ.-1) THEN
C           HAVE BEEN USING OLD PARTIALS, UPDATE THEM AND TRY AGAIN
            IWEVAL = MITER
            CFAIL = .TRUE.
            CALL F(N,T,Y,SAVE2,IPAR,RPAR,IERR)
            IF (IERR .NE. 0 ) GOTO 8000
            NFE = NFE + 1
            GO TO 110 
         ENDIF
         
         IJUS=0
         RED=0.5D+0
C    ***    failed at step 1 because of Newton
         GO TO 450
  
      ENDIF
  
      IWEVAL = -1
      HUSED = H
      NQUSED = NQ
      DO 140 I = 1,N
         Y(I,1) = (SAVE1(I)-ARH(I))
 140  CONTINUE
      DO 145 I=1,N
         SAVE2(I) = Y(I,1)*QQ
         Y(I,1) = SAVE1(I)
 145  CONTINUE
C     
C     UPDATE THE DIFFERENCES AT N+1
C
      DO 160 J = 2,L
         JM1 = J-1
         DO 150 I = 1,N
            Y(I,J) = Y(I,JM1) - YHOLD(I,JM1)
 150     CONTINUE
 160  CONTINUE
C     
C     COMPUTE ERROR IN THE SOLUTION
C
      DO 161 I=1,N
         AYI = DABS(Y(I,1))
         IF(ITOL.EQ.1) THEN
            SCALE(I) = YMAX(I)
         ELSE IF(ITOL.EQ.2) THEN
            SCALE(I) = RTOL(1)*AYI + ATOL(1)
         ELSE IF(ITOL.EQ.3) THEN
            SCALE(I) = RTOL(1)*AYI + ATOL(I)
         ELSE IF(ITOL.EQ.4) THEN
            SCALE(I) = RTOL(I)*AYI + ATOL(1)
         ELSE IF(ITOL.EQ.5) THEN
            SCALE(I) = RTOL(I)*AYI + ATOL(I)
         ENDIF
 161  CONTINUE
C ****
C ****  AMMEND
C ****  CHANGE 1,N BELOW TO 1,NVARS
C ****
      D = ZERO
      DO 170 I = 1,N
         D = D + ((Y(I,L)-YHOLD(I,L))/SCALE(I))**2
 170  CONTINUE
c     
c    STORING Y FROM FIRST STEP FOR USE IN THIRD STEP.
C
      IF(ITOL .EQ. 1) D = D/(RTOL(1)**2)
      IF(D.GT.E) GOTO 330
      DO 180 I = 1,N
         YNHOLD(I,1) = Y(I,1)
         YNHOLD(I,2) = SAVE2(I)
 180  CONTINUE
      
      KFAIL = 0
      IREDO = 0
      DO 190 I = 1,N
         ARH(I) = EL(2)*Y(I,1) 
 190  CONTINUE
      DO 210 J1 = 2,NQ
         JP1 = J1+1
         DO 200 I = 1,N
            ARH(I) = ARH(I) + EL(JP1)*Y(I,J1)
 200     CONTINUE
 210  CONTINUE
      
      CALL PRDICT(T,H,Y,L,N,SAVE2,NFE,F,IPAR,RPAR,IERR)
      IF (IERR .NE. 0 ) GOTO 8000
      DO 220 I = 1,N
         SAVE1(I) = Y(I,1)
         ERROR(I) = ZERO
 220  CONTINUE
      M2 = 0
C
C     FOR NOW WILL ASSUME THAT WE DO NOT WISH TO SAMPLE
C     AT THE N+2 STEP POINT
C      
      CALL ITRAT2(QQQ,Y,N,T,QI,BND,ARH,CRATE2,TCRAT2,M2,
     +     WORKED,YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,MF,
     +     1,ITOL,RTOL,ATOL,HUSED,NBSOL,NFE,NQUSED,F,
     +     RSP,IA,JA,P,IP,NSP,ISP,NONZ,LOUT,IPAR,RPAR,IERR)
      IF (IERR .NE. 0 ) GOTO 8000
      MEQC2 = MEQC2 + M2 + 1
C     
C       NOW CHECK TO SEE IF IT WAS SUCCESSFUL OR NOT
C
      IF (.NOT.WORKED) THEN
         NFAIL = NFAIL + 1
         IJUS=0
         RED=0.5D+0
C     ***have failed on step 2
         GOTO 450  
      ENDIF
C
C        IF WE ARE DOWN TO HERE THEN THINGS MUST HAVE CONVERGED
C
      LMP2=LMAX+2
      LMP3=LMAX+3
      DO 230 I = 1,N
         Y(I,LMP3) = (SAVE1(I)-ARH(I))
 230  CONTINUE
      DO 233 I=1,N
         Y(I,LMP2) = Y(I,LMP3)*QQ
         Y(I,LMP3) = SAVE1(I)
 233  CONTINUE
C
C     WE ARE NOW COMPUTING THE THIRD STAGE
C
      LL = L + 1
      T = TOLD + H
      DELST = ELST(1)-EL(1)
      NQP2 = NQ+2
      LMP2=LMAX+2
      DO 280 I=1,N
         SCALE(I) = H*(ELST(NQP2)*Y(I,LMP2)+DELST*YNHOLD(I,2))
         ARH(I) = SCALE(I)
         DO 270 J1 = 1,NQ
            ARH(I) = ARH(I) + ELST(J1+1)*YHOLD(I,J1)
 270     CONTINUE
 280  CONTINUE
      DO 290 I = 1,N
         SAVE2(I) = YNHOLD(I,2)
         Y(I,1) = YNHOLD(I,1)
 290  CONTINUE
      M3STEP = 0
 300  CONTINUE
      DO 310 I=1,N
         SAVE1(I) = (-Y(I,1)  + QI*SAVE2(I) + ARH(I))/qqq 
 310  CONTINUE
C
      IF (MF .EQ. 25 .OR. MF .EQ. 26) THEN 
         IFLAG = 0
         CALL CDRV(N,P,P,IP,IA,JA,PW,SAVE1,SAVE1,NSP,ISP,
     +        RSP,IEXTRA,4,IFLAG)
         IF (IFLAG .NE. 0) THEN
            WRITE(LOUT,4444) IFLAG
 4444       FORMAT('FALTAL ERROR IN CDRV ROUTINE 
     +           (SPARSE SOLVER) IFLAG RETURN IS ', I3)
            STOP            
         ENDIF
         NBSOL = NBSOL + 1
      ENDIF
      
      DO 321 I=1,N
         AYI = DABS(Y(I,1))
         IF(ITOL.EQ.1) THEN
            SCALE(I) = YMAX(I)
         ELSE IF(ITOL.EQ.2) THEN
            SCALE(I) = RTOL(1)*AYI + ATOL(1)
         ELSE IF(ITOL.EQ.3) THEN
            SCALE(I) = RTOL(1)*AYI + ATOL(I)
         ELSE IF(ITOL.EQ.4) THEN
            SCALE(I) = RTOL(I)*AYI + ATOL(1)
         ELSE IF(ITOL.EQ.5) THEN
            SCALE(I) = RTOL(I)*AYI + ATOL(I)
         ENDIF
 321  CONTINUE
      D = ZERO
      DO 320 I = 1,N
         D = D + (SAVE1(I)/SCALE(I))**2
         Y(I,1) = Y(I,1) + SAVE1(I)
 320  CONTINUE
      IF(ITOL .EQ. 1) D = D/(RTOL(1)**2)
      IF ((D*DMIN1(ONE,2.0D+0*CRATE1)).LE.BND) GO TO 360
      IF (M3STEP.EQ.4) THEN
c         WRITE (LOUT,9000)
         IJUS=1
         RED=0.5D+0
C     ****  step 3 fails
         NFAIL = NFAIL + 1
         GO TO 450
      ENDIF
      
      M3STEP = M3STEP + 1
      CALL F(N,T,Y,SAVE2,IPAR,RPAR,IERR)
      IF (IERR .NE. 0) GOTO 8000
      NFE = NFE + 1 
      GO TO 300
      
 330  KFAIL = KFAIL - 1
C **********************************************************************
C     THE ERROR TEST FAILED. KFAIL KEEPS TRACK OF MULTIPLE FAILURES.
C     RESTORE T AND THE Y ARRAY TO THEIR PREVIOUS VALUES AND PREPARE TO
C     TRY THE STEP AGAIN. COMPUTE THE OPTIMAL STEP SIZE FOR THIS ORDER
C     AND ONE ORDER LOWER.
C **********************************************************************
C     ***  failed on step 1 because of accuracy
C     COMPUTE ERROR IN THE SOLUTION
C
      NFAIL = NFAIL + 1
      DO 561 I=1,N
         AYI = DABS(Y(I,1))
         IF(ITOL.EQ.1) THEN
            SCALE(I) = YMAX(I)
         ELSE IF(ITOL.EQ.2) THEN
            SCALE(I) = RTOL(1)*AYI + ATOL(1)
         ELSE IF(ITOL.EQ.3) THEN
            SCALE(I) = RTOL(1)*AYI + ATOL(I)
         ELSE IF(ITOL.EQ.4) THEN
            SCALE(I) = RTOL(I)*AYI + ATOL(1)
         ELSE IF(ITOL.EQ.5) THEN
            SCALE(I) = RTOL(I)*AYI + ATOL(I)
         ENDIF
 561  CONTINUE
      DDOWN = ZERO
      TWODOWN = ZERO
      DO 1700 I=1,N
         DDOWN = DDOWN + ((Y(I,L))/SCALE(I))**2
         TWODOWN = TWODOWN + ((Y(I,l-1))/SCALE(I))**2 
 1700 CONTINUE
      IF(ITOL .EQ. 1) D = D/(RTOL(1)**2)
      T = TOLD
      HOLD = H
      IF(NQ.GT.1) FFAIL = 0.5D+0/DBLE(FLOAT(NQ))
      IF(NQ.GT.2) FRFAIL = 0.5D+0/DBLE(FLOAT(NQ-1))
      EFAIL = 0.5D+0/DBLE(FLOAT(L))
      CALL CPYARY(N*L,YHOLD,Y)
      RMAX = 2.0D+0 
      IF (DABS(H).LE.HMIN*1.00001D+0) THEN
C
C        REQUESTED ERROR NOT POSSIBLE WITH GIVEN HMIN
C
         KFLAG = -1 
         HOLD = H
         RETURN
  
      ENDIF
  
      IF (KFAIL.LE.-3) GO TO 340
      IREDO = 2
C
C     PREDICTING A NEW H AFTER INSUFFICIENT ACCURACY
C
      PRFAIL = ((D/(0.2D+0*E))**EFAIL)*1.5D+0 + 1.6D-6
      PLFAIL = ((DDOWN/(0.2D+0*EDN))**FFAIL)*1.5D+0+1.7D-6
      IF(NQ.GT.2) PLLFAIL =((TWODOWN/(0.2D+0*EDDN))**FRFAIL)*
     +     1.5D+0+1.7d-6
      IF(PLLFAIL.GT.PLFAIL) PLFAIL=PLLFAIL
      IF(PLFAIL.LT.PRFAIL.AND.NQ.NE.1) THEN
         NEWQ=NQ-1
         NQ=NEWQ
         RH=ONE/(PLFAIL*DBLE(FLOAT(-KFAIL)))
         L=NQ+1
         CALL COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD)
         RC=RC*EL(1)/OLDLO
         OLDLO=EL(1)
         CALL ERRORS(N,TQ,EDN,E,EUP,BND,EDDN)
      ELSE
         NEWQ = NQ
         RH = ONE/ (PRFAIL*DBLE(FLOAT(-KFAIL)))
      ENDIF
      GO TO 40
C **********************************************************************
C     CONTROL REACHES THIS STAGE IF 3 OR MORE FAILURES HAVE OCCURED.
C     IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE Y
C     ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST DERIVATIVE
C     IS RE-COMPUTED, AND THE ORDER IS SET TO 1. THEN H IS REDUCED BY A
C     FACTOR OF 10, AND THE STEP IS RETRIED. AFTER A TOTAL OF 7
C     FAILURES AN EXIT IS TAKEN WITH KFLAG=-2.
C **********************************************************************
 340  IF (KFAIL.EQ.-7) THEN
C        ERROR SMALLER THAN CAN BE HANDLED FOR PROBLEM
         KFLAG = -2 
         HOLD = H
         RETURN
  
      ENDIF

C     ********************************* 
C     START FROM ORDER 1 AGAIN    *
C     ********************************* 
      JCHANG = 1
      RH = DMAX1(HMIN/DABS(H),0.1D+0)
      CALL HCHOSE(RH,H,OVRIDE)
      H = H*RH
      CALL F(N,T,YHOLD,SAVE1,IPAR,RPAR,IERR) 
      IF (IERR .NE. 0 ) GOTO 8000
      NFE = NFE + 1 
      DO 350 I = 1,N
         Y(I,1) = YHOLD(I,1)
         Y(I,2) = H*SAVE1(I)
         YHOLD(I,2) = Y(I,2)
 350  CONTINUE
      IWEVAL = MITER
      CFAIL = .TRUE.
C     SINCE WE HAVE HAD PROBLEMS PROCEED WITH THIS ORDER
C     FOR 10 STEPS (IF WE CAN)
      IDOUB = 10
      IF (NQ.EQ.1) GO TO 60
      NQ = 1
      L = 2
C     RESET ORDER, RECALCULATE ERROR BOUNDS
      CALL COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD)
      LMAX = MAXDER + 1
      RC = RC*EL(1)/OLDLO
      OLDLO = EL(1) 
      CALL ERRORS(N,TQ,EDN,E,EUP,BND,EDDN)
C     NOW JUMP TO NORMAL CONTINUATION POINT
      GO TO 60
C **********************************************************************
C     THE ITERATION FOR THE CORRECTED SOLUTION HAS CONVERGED.
C     UPDATE THE Y ARRAY.
C **********************************************************************
 360  CONTINUE
C   ****
C   **** AMMEND ****
C   **** CHANGE 1,N BELOW TO 1,NVARS
C   ****
      DEMB=0.0D+0
      DO 361 I=1,N
         DEMB=DEMB+((Y(I,1)-YNHOLD(I,1))/SCALE(I))**2
 361  CONTINUE
      IF(DEMB.GT.4.0D+0*DBLE(FLOAT(N))) THEN
         IEMB=1
         IJUS=1
         RED=0.5D+0
C     ***  failed because of embedded error estimate
         NFAIL = NFAIL + 1
         GOTO 450
      ENDIF
      DO 380 J2 = 2,LL
         J2M1=J2-1
         DO 370 I = 1,N
            Y(I,J2) = Y(I,J2M1) - YHOLD(I,J2M1)
 370     CONTINUE
 380  CONTINUE
C ---------------------------------------------------------------------
C     IF THE COUNTER IDOUB EQUALS 2 AND WE ARE NOT ALREADY USING THE  
C     MAXIMUM ALLOWABLE ORDER , STORE Y(I,LMAX+4) WHICH IS USED IN
C     ASSESSING THE POSSIBILITY OF INCREASING THE ORDER. IF IDOUB = 0 
C     CONTROL PASSES TO 480 WHERE AN ATTEMPT TO CHANGE THE STEPSIZE AND
C     ORDER IS MADE.
C ----------------------------------------------------------------------
      IF (IDOUB.EQ.2.AND.L.NE.LMAX) THEN
         LMP4=LMAX+4
         DO 390 I = 1,N
            Y(I,LMP4) = Y(I,LL)
 390     CONTINUE
      ENDIF
      IDOUB = IDOUB - 1
      TRANGE=(TEND-TOLD-H)*H
      IF(TRANGE.LT.0.0D+0) THEN
         IDOUB = IDOUB + 2
         GOTO 440
      ENDIF
      JCHANG = 0
      IF (IDOUB.EQ.0) THEN
         SAMPLE = .FALSE.
         ISAMP = ISAMP + 1
         IF (ISAMP.EQ.4) THEN 
            SAMPLE = .TRUE.
            ISAMP = 0
         ENDIF
  
C **********************************************************************
C        NOW COMPUTE THE FACTORS PR1, PR2 AND PR3, BY WHICH 
C        H COULD BE DIVIDED AT ORDER NQ-1, ORDER NQ AND ORDER NQ+1
C        RESPECTIVELY. THE SMALLEST OF THESE IS DETERMINED AND THE NEW
C        ORDER CHOSEN ACCORDINGLY. IF THE ORDER IS TO BE INCREASED WE 
C        MUST COMPUTE ONE MORE BACKWARD DIFFERENCE.
C **********************************************************************
         PR3 = 1.D+20
         FAC = 1.5D+0
         IF(IEMB.EQ.1) FAC = 1.8D+0
         DO 400 I = 1,N
            AYI = DABS(Y(I,1))
            IF(ITOL.EQ.1) THEN
               VHOLD = YMAX(I)
            ELSE IF(ITOL.EQ.2) THEN 
               VHOLD = RTOL(1)*AYI + ATOL(1)
            ELSE IF(ITOL.EQ.3) THEN
               VHOLD = RTOL(1)*AYI + ATOL(I)
            ELSE IF(ITOL.EQ.4) THEN
               VHOLD = RTOL(I)*AYI + ATOL(1)
            ELSE IF(ITOL.EQ.5) THEN
               VHOLD = RTOL(I)*AYI + ATOL(I)
            ENDIF
            SCALE(I)=VHOLD
 400     CONTINUE
         IF(L.NE.LMAX) THEN
            LMP4 = LMAX + 4
            DUP = ZERO
            DO 401 I=1,N
               DUP = DUP + ((Y(I,LL)-Y(I,LMP4))/SCALE(I))**2
 401        CONTINUE
            IF(ITOL .EQ. 1) DUP = DUP/(RTOL(1)**2)
            ENQ3 = 0.5D+0/DBLE(FLOAT(L+1))
            PR3 = ((DUP/EUP)**ENQ3)*(FAC+0.2D+0) + 1.8D-6
         ENDIF
         
         ENQ2 = 0.5D+0/DBLE(FLOAT(L))
         D = ZERO
         DDOWN=ZERO
         DO 410 I = 1,N
            D = D + (Y(I,LL)/SCALE(I))**2
 410     CONTINUE
         IF(ITOL .EQ.1) D = D/(RTOL(1)**2)
         PR2 = ((D/E)**ENQ2)*FAC + 1.6D-6
         PR1 = 1.D+20
         IF (NQ.GT.1) THEN
            DDDOWN=ZERO
            DDOWN = ZERO
            DO 57420 I = 1,N
               DDOWN = DDOWN + (Y(I,L)/SCALE(I))**2 
               DDDOWN = DDDOWN + (Y(I,L-1)/SCALE(I))**2
57420       CONTINUE
            IF(ITOL .EQ. 1) DDOWN = DDOWN/(RTOL(1)**2)
            ENQ1 = 0.5D+0/DBLE(FLOAT(NQ))
            PR1 = ((DDOWN/EDN)**ENQ1)*(FAC+0.1D+0) + 1.7D-6
            IF(NQ.GT.2) THEN
               ENQ0 = 0.5D+0/DBLE(FLOAT(NQ-1))
               PR0 = ((DDDOWN/EDDN)**ENQ0)*(FAC+0.1D+0) + 1.7D-6
               IF(PR0.GT.PR1) PR1 = PR0
               IF(DDDOWN.LT.DDOWN) DDOWN = DDDOWN
            ENDIF
         ENDIF
         IF(L.EQ.LMAX) DUP = 0.0D+0
         IF(NQ.LE.1) GOTO 6578
         IF(DUP.GT.D.AND.D.GT.DDOWN) THEN
            PR2=1.0D+30
            PR3=1.0D+30
         ENDIF
 6578    CONTINUE
         IF (PR2.LE.PR3) THEN 
            IF (PR2.GT.PR1) THEN
               NEWQ = NQ - 1
               RH = 1.0D+0/PR1  
            ELSE
               NEWQ = NQ
               RH = 1.0D+0/PR2
            ENDIF
         ELSE IF (PR3.LT.PR1) THEN
            NEWQ = L
            RH = 1.0D+0/PR3
         ELSE
            NEWQ = NQ - 1
            RH = 1.0D+0/PR1
         ENDIF
         IEMB=0
         IF(RH.GT.1.0D+0.AND.RH.LT.1.1D+0) THEN
            IDOUB=10
            NQ=NQUSED
            L=NQ+1
            GOTO 440
         ENDIF
         RH = DMIN1(RH,RMAX)
         CALL HCHOSE(RH,H,OVRIDE)
         IF ((JSINUP.LE.20).AND.(KFLAG.EQ.0).AND.(RH.LT.1.1D+0)) THEN 
C           WE HAVE RUN INTO PROBLEMS
            IDOUB = 10
            NQ = NQUSED
            L = NQ + 1
            GO TO 440
            
         ENDIF
C **********************************************************************
C        IF THERE IS A CHANGE IN ORDER, RESET NQ, L AND THE 
C        COEFFICIENTS. IN ANY CASE H IS RESET  AND THE
C        Y ARRAY IS RE-SCALED 
C **********************************************************************
         IF (NEWQ.NE.NQ) THEN 
            IF (NEWQ.GT.NQ) THEN
C              ADD AN EXTRA TERM TO THE HISTORY ARRAY
               DO 430 I = 1,N 
                  Y(I,LL) = Y(I,L) - YHOLD(I,L)
 430           CONTINUE
            ENDIF
  
            NQ = NEWQ
            L = NQ + 1
C     RESET ORDER,RECALCULATE ERROR BOUNDS
            CALL COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD)
            LMAX = MAXDER + 1 
            RC = RC*EL(1)/OLDLO
            OLDLO = EL(1)
            CALL ERRORS(N,TQ,EDN,E,EUP,BND,EDDN)
         ENDIF
  
         RH = DMAX1(RH,HMIN/DABS(H))
         RH = DMIN1(RH,HMAX/DABS(H),RMAX)
         CALL RSCALE(N,L,RH,Y)
         RMAX = 10.0D+0
         JCHANG = 1 
         H = H*RH
         RC = RC*RH 
         IF(JSNOLD.GT.IBND) RC=ZERO
         IDOUB = L + 1
      ENDIF
      
  440 CONTINUE
C ----------------------------------------------------------------------
C     STORE THE Y ARRAY IN THE MATRIX YHOLD.  STORE IN THE Y ARRAY THE
C     INFORMATION NECESSARY TO PERFORM AN INTERPOLATION TO FIND THE
C     SOLUTION AT THE SPECIFIED OUTPUT POINT IF APPROPRIATE.
C ----------------------------------------------------------------------
      CALL CPYARY(N*L,Y,YHOLD)
      NSTEP = NSTEP + 1
      JSINUP = JSINUP + 1
      JSNOLD = JSNOLD + 1
      JSTART = NQUSED
      T = TOLD + HUSED
      HOLD = H
      KFAIL = 0
      NEWPAR = 0
      CFAIL = .FALSE.
      RETURN
 450  CONTINUE
      FINISH = .FALSE.
      T=TOLD
      RMAX=2.0D+0
      DO 460 J1=1,L 
         DO 460 I=1,N
            Y(I,J1)=YHOLD(I,J1) 
 460     CONTINUE
      IF(DABS(H).LE.HMIN*1.00001D+0) THEN
C     
C   CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED
C
         IF(NSTEP.EQ.0) THEN 
            KFLAG=-1
         ELSE
            KFLAG=-3
         ENDIF
C
C    TO SUPPRESS ERROR MESSAGES AT START AS H MAY 
C    HAVE BEEN TOO LARGE ON THE FIRST STEP.
C
         HOLD=H
         FINISH = .TRUE.
      ENDIF
      RH = RED
      IREDO=1
C
C     TRY AGAIN WITH UPDATED PARTIALS
C
 8000 IF (IERR .NE. 0 ) RETURN
      IF(IJUS.EQ.0) CALL HCHOSE(RH,H,OVRIDE)
      IF(.NOT.FINISH) THEN
         GO TO 40
      ELSE
         RETURN
      ENDIF
 9000 FORMAT (1X,' CORRECTOR HAS NOT CONVERGED')
      END 
C ------------------- END OF SUBROUTINE STIFF --------------------------

      SUBROUTINE RSCALE(N,L,RH,Y)
      
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C     .. SCALAR ARGUMENTS ..
      INTEGER L,N
C     ..
C     .. ARRAY ARGUMENTS ..
      DIMENSION  Y(N,12)
C     ..
C     .. LOCAL SCALARS ..
      INTEGER I,J,J1
C     ..
C     .. LOCAL ARRAYS ..
      DIMENSION  DI(8,8)
C     ..
C     .. DATA STATEMENTS ..
  
C     SUBROUTINE IS FOR RESCALING THE HISTORY ARRAY AFTER A CHANGE IN 
C     STEPSIZE
C
C     N      ORDER OF THE PROBLEM
C     L      NUMBER OF TERMS IN THE HISTORY ARRAY TO BE RESCALED
C     RH     RATIO OF THE STEPSIZE CHANGE (I.E. RH = HNEW/HOLD)
C     Y()    THE HISTORY ARRAY
C
      DATA  ZERO/0.0D+0/
C     ..
 
      DI(2,2) = RH
      IF (L.GT.2) THEN
         TA = RH*RH 
         DI(2,3) = RH* (1.0D+0-RH)/2.0D+0
         DI(3,3) = TA
         IF (L.GT.3) THEN
            TB = TA*RH
            DI(2,4) = RH* ((RH-3.0D+0)*RH+2.0D+0)/6.0D+0
            DI(3,4) = TA* (1.0D+0-RH)
            DI(4,4) = TB
            IF (L.GT.4) THEN
               TC = TB*RH
               DI(2,5) = - (((RH-6.0D+0)*RH+11.0D+0)*RH-6.0D+0)*RH/
     +                      24.0D+0
               DI(3,5) = TA* ((7.0D+0*RH-18.0D+0)*RH+11.0D+0)/12.0D+0 
               DI(4,5) = 1.5D+0*TB* (1.0D+0-RH)
               DI(5,5) = TC
               IF (L.GT.5) THEN
                  TD = TC*RH
                  DI(2,6) = ((((RH-10.0D+0)*RH+35.0D+0)*RH-50.0D+0)
     +                        *RH+24.0D+0)*RH/120.0D+0
                  DI(3,6) = - (((3.0D+0*RH-14.0D+0)*RH+21.0D+0)*RH
     +                       -10.0D+0)*TA/12.0D+0 
                  DI(4,6) = ((5.0D+0*RH-12.0D+0)*RH+7.0D+0)*TB/4.0D+0 
                  DI(5,6) = 2.0D+0*TC* (1.0D+0-RH)
                  DI(6,6) = TD
                  IF (L.GT.6) THEN
                     TE = TD*RH
                     DI(2,7) = -RH* (RH-1.0D+0)* (RH-2.0D+0)*
     +                          (RH-3.0D+0)*(RH-4.0D+0)*(RH-5.0D+0)/
     +                           720.0D+0
                     DI(3,7) = TA* ((((62.0D+0*RH-450.0D+0)*RH+
     +                         1190.0D+0)*RH-1350.0D+0)*RH+548.0D+0)
     +                         /720.0D+0
                     DI(4,7) = TB* (((-18.0D+0*RH+75.0D+0)*RH
     +                          -102.0D+0)*RH+45.0D+0)/24.0D+0
                     DI(5,7) = TC* ((13.0D+0*RH-30.0D+0)*RH+17.0D+0)
     +                         /6.0D+0
                     DI(6,7) = 2.5D+0*TD* (1.0D+0-RH)
                     DI(7,7) = TE
                     IF (L.GT.7) THEN
                        TF = TE*RH
                        DI(2,8) = RH*(RH-1.0D+0)*(RH-2.0D+0)*(RH
     +                            -3.0D+0)*(RH-4.0D+0)*(RH-5.0D+0)
     +                            *(RH-6.0D+0)/5040.0D+0
                        DI(3,8) = TA* ((((((-126.0D+0*RH)+1302.0D+0)*RH-
     +                            5250.0D+0)*RH+10290.0D+0)*RH-9744.0D+0
     +                            )*RH+3528.0D+0)/5040.0D+0 
                        DI(4,8) = TB* ((((43.0D+0*RH-270.0D+0)*RH+
     +                            625.0D+0)*RH-630.0D+0)*RH+232.0D+0) 
     +                           /120.0D+0
                        DI(5,8) = TC* (((-10.0D+0*RH+39.0D+0)*RH-
     +                            50.0D+0)*RH+21.0D+0)/6.0D+0
                        DI(6,8) = TD* ((20.0D+0*RH-45.0D+0)*RH+25.0D+0
     +                            )/6.0D+0
                        DI(7,8) = 3.0D+0*TE* (1.0D+0-RH)
                        DI(8,8) = TF
                     ENDIF
  
                  ENDIF
  
               ENDIF
  
            ENDIF
  
         ENDIF
  
      ENDIF
  
      DO 30 I = 1,N 
        DO 20 J = 2,L
          ZZ = ZERO 
          DO 10 J1 = J,L
            ZZ = ZZ + DI(J,J1)*Y(I,J1)
   10     CONTINUE
          Y(I,J) = ZZ
   20   CONTINUE
   30 CONTINUE
      RETURN
  
      END 
C-----------------------------------------------------------------------------

      SUBROUTINE CPYARY(NELEM,SOURCE,TARGET)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     COPIES THE ARRAY SOURCE() INTO THE ARRAY TARGET()
C
C     THIS SUBROUTINE COULD BE REPLACED BY THE BLAS ROUTINE SCOPY
C     (AFTER CHANGING THE ARGUMENT LIST APPROPRIATELY)
C
C     .. SCALAR ARGUMENTS ..
      INTEGER NELEM 
C     ..
C     .. ARRAY ARGUMENTS ..
      DIMENSION  SOURCE(NELEM),TARGET(NELEM)
C     ..
C     .. LOCAL SCALARS ..
      INTEGER I
C     ..
      DO 10 I = 1,NELEM
         TARGET(I) = SOURCE(I) 
 10   CONTINUE
      RETURN      
      END 

C-----------------------------------------------------------------------------

      SUBROUTINE HCHOSE(RH,H,OVRIDE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)

      COMMON / STPSZE / HSTPSZ(2,14)
      LOGICAL OVRIDE
C     
C     FIRST MOVE ALL ELEMENTS DOWN ONE PLACE
C
      IF (H.NE.HSTPSZ(2,1)) THEN
         DO 5 I=12,2,-1
           I2=I-1
           HSTPSZ(1,I)=HSTPSZ(1,I2)
    5      HSTPSZ(2,I)=HSTPSZ(2,I2)
C
C          NOW INSERT VALUE OF H USED BEFORE THIS CALL
C
         HSTPSZ(1,2)=H/HSTPSZ(2,1)
         HSTPSZ(2,1)=H
      ENDIF
C
C     NOW DECIDE ON THE NEW CHANGE
C
      IF (RH.GT.1.0D+0) THEN
         OVRIDE=.FALSE.
      ELSE IF (HSTPSZ(1,2).LE.1.0D+0) THEN
         OVRIDE=.FALSE.
      ELSE IF ((RH*H).LE.HSTPSZ(2,2)) THEN
         OVRIDE=.FALSE.
      ELSE
         RH=HSTPSZ(2,2)/H 
         OVRIDE=.TRUE.
      ENDIF
      HSTPSZ(1,1)=RH
      RETURN
      END 
C-----------------------------------------------------------------------------
C     THE NEXT ROUTINE PREFORMS PREPROCESSING RELATED TO THE SPARSE LINEAR
C     SYSTEMS.
C     THE OPERATIONS THAT ARE PERFORMED HERE ARE  
C  * COMPUTE SPARSENESS STRUCTURE OF THE JACOBIAN MATRIX ACCORDING TO MF,
C  * COMPUTE GROUPING OF COLUMN INDICES (MF = 26),
C  * COMPUTE A NEW ORDERING OF ROWS AND COLUMNS OF THE MATRIX,
C  * REORDER JA CORRESPONDING TO THE NEW ORDERING,
C  * PERFORM A SYMBOLIC LU FACTORIZATION OF THE MATRIX.
C-----------------------------------------------------------------------------
C
      SUBROUTINE SPARSDET(N,T,Y0,PWCOPY,SCALE,IA,JA,SAVF,P,IP,NSP,
     +     ISP,RSP,NONZ,NGP,JGP,IGP,MF,INCL,JDONE,LOUT,F,PDERV,NNZ,
     +     IPAR,RPAR,IERR)
C 
C     ..ARRAY ARGUMENTS ..
      DOUBLE PRECISION T
      DOUBLE PRECISION PWCOPY(NONZ),Y0(N),RSP(NSP),SAVF(N),YJ
      DOUBLE PRECISION FAC,SCALE(N),FTEM(1),DQ,RPAR(*)
      INTEGER P(N),IP(N),IA(N+1),JA(NONZ),ISP(NSP), JGP(N),IGP(NONZ),
     +     INCL(N),JDONE(N),IPAR(*)
C     
      LOGICAL ORDER
C     .. EXTERNAL SUBROUTINES ..
      EXTERNAL F,PDERV,ODRV,CDRV
C     ..
C     .. INTRINSIC FUNCTIONS ..
      INTRINSIC DABS,DMAX1
C     ..
C     .. COMMON BLOCKS ..  
      COMMON /ORDERING/ORDER      
C...
C...  PERTURB THE INITIAL CONDITIONS, Y0, FOR STRUCTURE DETERMINATION.      
C...      
      DO I=1,N             
         FAC= 1.0D0+1.0D0/(DFLOAT(I)+1.D0)
         Y0(I) =Y0(I) +FAC*DSIGN(SCALE(I),Y0(I))          
      ENDDO
C...  
C...  COMPUTE THE SPARSITY STRUCTURE FROM THE USER SUPPLIES JACOBIAN 
C...  ROUTINE PDERV.      
C...  
      IF (MF .EQ. 25) THEN
         KK = 1
         IA(1) =1         
         DO 111 J =1,N
            JA(KK) = J
            KK = KK+ 1                        
            DO K=1,N
               SAVF(K) = 0.0D0
            ENDDO             
            CALL PDERV(N,T,Y0,J,SAVF,IPAR,RPAR,IERR)            
            DO 112 II=1,N   
               IF (DABS(SAVF(II)) .LE. 0.0D0 ) GOTO 112               
               IF (II .EQ. J) GOTO 112               
               JA(KK) = II               
               KK=KK+1
 112        CONTINUE
            IA(1+J) = KK
 111     CONTINUE
      ELSE
C....  
C.... MF = 26  COMPUTE SPARSITY STRUCTURE FROM RESULTS OF N+1 CALLS TO F
C....
 100     KK = 1
         IA(1) = 1      
         CALL F(N,T,Y0,SAVF,IPAR,RPAR,IERR)
         DO 120 J = 1,N
            JA(KK) = J
            KK = KK + 1
            YJ = Y0(J)
            DYJ = DSIGN(SCALE(J),YJ)
            Y0(J) = YJ + DYJ
            CALL F(N,T,Y0,PWCOPY,IPAR,RPAR,IERR)
            Y0(J) = YJ            
            DO 110 II = 1,N
               DQ = (PWCOPY(II) - SAVF(II))/DYJ               
               IF (DABS(DQ) .LE. 0.D0) GOTO 110
               IF (II .EQ. J) GOTO 110
               JA(KK) = II
               KK = KK + 1
 110        CONTINUE
            IA(1+J) = KK
 120     CONTINUE
      ENDIF
C.....
      NNZ = IA(1+N)-1
      IF (NNZ .GT. NONZ) THEN 
         WRITE(LOUT,*) 'NUMBER OF NONZEROS ',NNZ, 'NOT',NONZ
         STOP
      ENDIF
C....
C.... COMPUTE COLUMN GROUPING FOR SPARSE NUMERICAL JACOBIAN I.E. MF = 26
C....                                
      IF ( MF .EQ. 26) THEN
         MAXG = N+1  
         CALL JGROUP(N,IA,JA,MAXG,NGP,IGP,JGP,INCL,JDONE,IER,NNZ)
      ENDIF
C.... 
C.... COMPUTE NEW ORDERING OF NONZERO  ENTERIS IN THE ITERATION MATRIX
C....
      DO I=1,N
         P(I) = I               
      ENDDO
      CALL ODRV(N,IA,JA,PWCOPY,P,IP,NSP,ISP,1,IFLAG)
      IF (IFLAG .NE. 0) THEN  
         WRITE(LOUT,*)'FATAL ERROR IN THE ODRV ROUTINE (SPARSE ORDERING)
     +        IFALG RETURN IS= ',IFLAG
         STOP            
      ENDIF          
C.... 
C.... REORDER AND DO SYMBOLIC LU FACTORIZATION OF THE ITERATION MATRIX
C.... 
      DO I=1, NNZ
         PWCOPY(I) =0.0D0               
      ENDDO
      CALL CDRV(N,P,P,IP,IA,JA,PWCOPY,FTEM,FTEM,NSP,ISP,
     +     RSP,IEXTRA,5,IFLAG)
      IF (IFLAG .EQ. 0) THEN 
         ORDER = .TRUE.
      ELSE                  
         WRITE(LOUT,1111) IFLAG
 1111    FORMAT('FATAL ERROR IN THE CDRV ROUTINE (SYMBOLIC 
     +        FACTORIZATION) IFLAG RETURN IS ', I3)
         STOP              
      ENDIF
C
      RETURN
      END
C-----------------------------------------------------------------------------

      SUBROUTINE JGROUP(N,IA,JA,MAXG,NGRP,IGP,JGP,INCL,JDONE,IER,NNZ)      
      INTEGER N, IA, JA, MAXG, NGRP, IGP, JGP, INCL, JDONE, IER, NNZ 
      DIMENSION IA(1), JA(1), IGP(1), JGP(N), INCL(N), JDONE(N)
C-----------------------------------------------------------------------------
C THIS SUNROUTINE CONSTRUCTS GROUPING OF THE COLUMN INDICES OF THE
C JACOBIAN MATRIX, USED IN THE NUMERICAL EVALUATION OF THE JACOBIAN BY
C FINITE DIFFERENCES.     
C 
C INPUT.. 
C N      = THE ORDER OF THE MATRIX. 
C IA,JA  = SPARSE STRUCURE DESCRIPTORS OF THE MATRIX BY ROWS.
C MAXG   = LENGTH OF AVAILABLE STORAGE IN THE IGP ARRAY.
C
C OUTPUT..
C NGRP   = NUMBER OF GROUPS.
C JGP    = ARRAY OF LENGTH N CONTAINING THE COLUMN INDICES BY GROUPS.
C IGP    = POINTER ARRAY OF LENGTH NGRP + 1 TO THE LOCATIONS IN JGP OF
C          THE BEGINNING OF EACH GROUP.
C IER    = ERROR INDICATOR. IER = 0 IF NO ERROR OCCURRED, OR 1 IF
C          MAXG WAS INSUFFICIENT.
C
C INCL AND JDONE ARE WORKING ARRAYS OF LENGTH N.
C-----------------------------------------------------------------------------
      INTEGER I, J, K, KMIN, KMAX, NCOL, NG 
C
      IER = 0
      DO 10 JK = 1,N
 10      JDONE(JK) = 0
      NCOL = 1
      DO 60 NG = 1,MAXG
         IGP(NG) = NCOL
         DO 20 I = 1,N
 20         INCL(I) = 0
         DO 50 J = 1,N
C...
C...  REJECT COLUMN J IF IT IS ALREADY IN A GROUP.
C...
            IF (JDONE(J) .EQ. 1) GOTO 50
            KMIN = IA(J)
            KMAX = IA(J+1) - 1
            DO 30 K = KMIN, KMAX
C...
C...  REJECT COLUMN J IF IT OVERLAPS ANY COLUMN ALREADY IN THIS GROUP.
C...
               I = JA(K)           
               IF (INCL(I) .EQ. 1) GOTO 50
 30         CONTINUE
C...
C...  ACCEPT COLUMN J INTO GROUP NG.
C...
            JGP(NCOL) = J          
            NCOL = NCOL + 1          
            JDONE(J) = 1           
            DO 40 K = KMIN,KMAX
               I = JA(K)
               INCL(I) = 1
 40         CONTINUE
 50      CONTINUE          
C...
C...  STOP IF THIS GROUP IS EMPTY (GROUPING IS COMPLETE).
C...
         IF (NCOL .EQ. IGP(NG)) GOTO 70
 60   CONTINUE
C...
C...  ERROR RETURN IF NOT ALL COLUMNS WERE CHOSEN (MAXG TOO SMALL).
C...
      IF (NCOL.LE. N) GOTO 80
      NG = MAXG
 70   NGRP = NG - 1               
      RETURN
 80   IER = 1
      RETURN
      END
C----------------------- END OF SUBROUTINE JGROUP  ----------------------

      DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. Scalar Arguments ..
      CHARACTER          CMACH
*     ..
*
*  Purpose
*  =======
*
*  DLAMCH determines double precision machine parameters.
*
*  Arguments
*  =========
*
*  CMACH   (input) CHARACTER*1
*          Specifies the value to be returned by DLAMCH:
*          = 'E' or 'e',   DLAMCH := eps
*          = 'S' or 's ,   DLAMCH := sfmin
*          = 'B' or 'b',   DLAMCH := base
*          = 'P' or 'p',   DLAMCH := eps*base
*          = 'N' or 'n',   DLAMCH := t
*          = 'R' or 'r',   DLAMCH := rnd
*          = 'M' or 'm',   DLAMCH := emin
*          = 'U' or 'u',   DLAMCH := rmin
*          = 'L' or 'l',   DLAMCH := emax
*          = 'O' or 'o',   DLAMCH := rmax
*
*          where
*
*          eps   = relative machine precision
*          sfmin = safe minimum, such that 1/sfmin does not overflow
*          base  = base of the machine
*          prec  = eps*base
*          t     = number of (base) digits in the mantissa
*          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise
*          emin  = minimum exponent before (gradual) underflow
*          rmin  = underflow threshold - base**(emin-1)
*          emax  = largest exponent before overflow
*          rmax  = overflow threshold  - (base**emax)*(1-eps)
*
* =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE, ZERO
      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            FIRST, LRND
      INTEGER            BETA, IMAX, IMIN, IT
      DOUBLE PRECISION   BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
     $                   RND, SFMIN, SMALL, T
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     ..
*     .. External Subroutines ..
      EXTERNAL           DLAMC2
*     ..
*     .. Save statement ..
      SAVE               FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN,
     $                   EMAX, RMAX, PREC
*     ..
*     .. Data statements ..
      DATA               FIRST / .TRUE. /
*     ..
*     .. Executable Statements ..
*
      IF( FIRST ) THEN
         FIRST = .FALSE.
         CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )
         BASE = BETA
         T = IT
         IF( LRND ) THEN
            RND = ONE
            EPS = ( BASE**( 1-IT ) ) / 2
         ELSE
            RND = ZERO
            EPS = BASE**( 1-IT )
         END IF
         PREC = EPS*BASE
         EMIN = IMIN
         EMAX = IMAX
         SFMIN = RMIN
         SMALL = ONE / RMAX
         IF( SMALL.GE.SFMIN ) THEN
*
*           Use SMALL plus a bit, to avoid the possibility of rounding
*           causing overflow when computing  1/sfmin.
*
            SFMIN = SMALL*( ONE+EPS )
         END IF
      END IF
*
      IF( LSAME( CMACH, 'E' ) ) THEN
         RMACH = EPS
      ELSE IF( LSAME( CMACH, 'S' ) ) THEN
         RMACH = SFMIN
      ELSE IF( LSAME( CMACH, 'B' ) ) THEN
         RMACH = BASE
      ELSE IF( LSAME( CMACH, 'P' ) ) THEN
         RMACH = PREC
      ELSE IF( LSAME( CMACH, 'N' ) ) THEN
         RMACH = T
      ELSE IF( LSAME( CMACH, 'R' ) ) THEN
         RMACH = RND
      ELSE IF( LSAME( CMACH, 'M' ) ) THEN
         RMACH = EMIN
      ELSE IF( LSAME( CMACH, 'U' ) ) THEN
         RMACH = RMIN
      ELSE IF( LSAME( CMACH, 'L' ) ) THEN
         RMACH = EMAX
      ELSE IF( LSAME( CMACH, 'O' ) ) THEN
         RMACH = RMAX
      END IF
*
      DLAMCH = RMACH
      RETURN
*
*     End of DLAMCH
*
      END
*
************************************************************************
*
      SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. Scalar Arguments ..
      LOGICAL            IEEE1, RND
      INTEGER            BETA, T
*     ..
*
*  Purpose
*  =======
*
*  DLAMC1 determines the machine parameters given by BETA, T, RND, and
*  IEEE1.
*
*  Arguments
*  =========
*
*  BETA    (output) INTEGER
*          The base of the machine.
*
*  T       (output) INTEGER
*          The number of ( BETA ) digits in the mantissa.
*
*  RND     (output) LOGICAL
*          Specifies whether proper rounding  ( RND = .TRUE. )  or
*          chopping  ( RND = .FALSE. )  occurs in addition. This may not
*          be a reliable guide to the way in which the machine performs
*          its arithmetic.
*
*  IEEE1   (output) LOGICAL
*          Specifies whether rounding appears to be done in the IEEE
*          'round to nearest' style.
*
*  Further Details
*  ===============
*
*  The routine is based on the routine  ENVRON  by Malcolm and
*  incorporates suggestions by Gentleman and Marovich. See
*
*     Malcolm M. A. (1972) Algorithms to reveal properties of
*        floating-point arithmetic. Comms. of the ACM, 15, 949-951.
*
*     Gentleman W. M. and Marovich S. B. (1974) More on algorithms
*        that reveal properties of floating point arithmetic units.
*        Comms. of the ACM, 17, 276-277.
*
* =====================================================================
*
*     .. Local Scalars ..
      LOGICAL            FIRST, LIEEE1, LRND
      INTEGER            LBETA, LT
      DOUBLE PRECISION   A, B, C, F, ONE, QTR, SAVEC, T1, T2
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMC3
      EXTERNAL           DLAMC3
*     ..
*     .. Save statement ..
      SAVE               FIRST, LIEEE1, LBETA, LRND, LT
*     ..
*     .. Data statements ..
      DATA               FIRST / .TRUE. /
*     ..
*     .. Executable Statements ..
*
      IF( FIRST ) THEN
         FIRST = .FALSE.
         ONE = 1
*
*        LBETA,  LIEEE1,  LT and  LRND  are the  local values  of  BETA,
*        IEEE1, T and RND.
*
*        Throughout this routine  we use the function  DLAMC3  to ensure
*        that relevant values are  stored and not held in registers,  or
*        are not affected by optimizers.
*
*        Compute  a = 2.0**m  with the  smallest positive integer m such
*        that
*
*           fl( a + 1.0 ) = a.
*
         A = 1
         C = 1
*
*+       WHILE( C.EQ.ONE )LOOP
   10    CONTINUE
         IF( C.EQ.ONE ) THEN
            A = 2*A
            C = DLAMC3( A, ONE )
            C = DLAMC3( C, -A )
            GO TO 10
         END IF
*+       END WHILE
*
*        Now compute  b = 2.0**m  with the smallest positive integer m
*        such that
*
*           fl( a + b ) .gt. a.
*
         B = 1
         C = DLAMC3( A, B )
*
*+       WHILE( C.EQ.A )LOOP
   20    CONTINUE
         IF( C.EQ.A ) THEN
            B = 2*B
            C = DLAMC3( A, B )
            GO TO 20
         END IF
*+       END WHILE
*
*        Now compute the base.  a and c  are neighbouring floating point
*        numbers  in the  interval  ( beta**t, beta**( t + 1 ) )  and so
*        their difference is beta. Adding 0.25 to c is to ensure that it
*        is truncated to beta and not ( beta - 1 ).
*
         QTR = ONE / 4
         SAVEC = C
         C = DLAMC3( C, -A )
         LBETA = C + QTR
*
*        Now determine whether rounding or chopping occurs,  by adding a
*        bit  less  than  beta/2  and a  bit  more  than  beta/2  to  a.
*
         B = LBETA
         F = DLAMC3( B / 2, -B / 100 )
         C = DLAMC3( F, A )
         IF( C.EQ.A ) THEN
            LRND = .TRUE.
         ELSE
            LRND = .FALSE.
         END IF
         F = DLAMC3( B / 2, B / 100 )
         C = DLAMC3( F, A )
         IF( ( LRND ) .AND. ( C.EQ.A ) )
     $      LRND = .FALSE.
*
*        Try and decide whether rounding is done in the  IEEE  'round to
*        nearest' style. B/2 is half a unit in the last place of the two
*        numbers A and SAVEC. Furthermore, A is even, i.e. has last  bit
*        zero, and SAVEC is odd. Thus adding B/2 to A should not  change
*        A, but adding B/2 to SAVEC should change SAVEC.
*
         T1 = DLAMC3( B / 2, A )
         T2 = DLAMC3( B / 2, SAVEC )
         LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
*
*        Now find  the  mantissa, t.  It should  be the  integer part of
*        log to the base beta of a,  however it is safer to determine  t
*        by powering.  So we find t as the smallest positive integer for
*        which
*
*           fl( beta**t + 1.0 ) = 1.0.
*
         LT = 0
         A = 1
         C = 1
*
*+       WHILE( C.EQ.ONE )LOOP
   30    CONTINUE
         IF( C.EQ.ONE ) THEN
            LT = LT + 1
            A = A*LBETA
            C = DLAMC3( A, ONE )
            C = DLAMC3( C, -A )
            GO TO 30
         END IF
*+       END WHILE
*
      END IF
*
      BETA = LBETA
      T = LT
      RND = LRND
      IEEE1 = LIEEE1
      RETURN
*
*     End of DLAMC1
*
      END
*
************************************************************************
*
      SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. Scalar Arguments ..
      LOGICAL            RND
      INTEGER            BETA, EMAX, EMIN, T
      DOUBLE PRECISION   EPS, RMAX, RMIN
*     ..
*
*  Purpose
*  =======
*
*  DLAMC2 determines the machine parameters specified in its argument
*  list.
*
*  Arguments
*  =========
*
*  BETA    (output) INTEGER
*          The base of the machine.
*
*  T       (output) INTEGER
*          The number of ( BETA ) digits in the mantissa.
*
*  RND     (output) LOGICAL
*          Specifies whether proper rounding  ( RND = .TRUE. )  or
*          chopping  ( RND = .FALSE. )  occurs in addition. This may not
*          be a reliable guide to the way in which the machine performs
*          its arithmetic.
*
*  EPS     (output) DOUBLE PRECISION
*          The smallest positive number such that
*
*             fl( 1.0 - EPS ) .LT. 1.0,
*
*          where fl denotes the computed value.
*
*  EMIN    (output) INTEGER
*          The minimum exponent before (gradual) underflow occurs.
*
*  RMIN    (output) DOUBLE PRECISION
*          The smallest normalized number for the machine, given by
*          BASE**( EMIN - 1 ), where  BASE  is the floating point value
*          of BETA.
*
*  EMAX    (output) INTEGER
*          The maximum exponent before overflow occurs.
*
*  RMAX    (output) DOUBLE PRECISION
*          The largest positive number for the machine, given by
*          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point
*          value of BETA.
*
*  Further Details
*  ===============
*
*  The computation of  EPS  is based on a routine PARANOIA by
*  W. Kahan of the University of California at Berkeley.
*
* =====================================================================
*
*     .. Local Scalars ..
      LOGICAL            FIRST, IEEE, IWARN, LIEEE1, LRND
      INTEGER            GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
     $                   NGNMIN, NGPMIN
      DOUBLE PRECISION   A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
     $                   SIXTH, SMALL, THIRD, TWO, ZERO
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMC3
      EXTERNAL           DLAMC3
*     ..
*     .. External Subroutines ..
      EXTERNAL           DLAMC1, DLAMC4, DLAMC5
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX, MIN
*     ..
*     .. Save statement ..
      SAVE               FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
     $                   LRMIN, LT
*     ..
*     .. Data statements ..
      DATA               FIRST / .TRUE. / , IWARN / .FALSE. /
*     ..
*     .. Executable Statements ..
*
      IF( FIRST ) THEN
         FIRST = .FALSE.
         ZERO = 0
         ONE = 1
         TWO = 2
*
*        LBETA, LT, LRND, LEPS, LEMIN and LRMIN  are the local values of
*        BETA, T, RND, EPS, EMIN and RMIN.
*
*        Throughout this routine  we use the function  DLAMC3  to ensure
*        that relevant values are stored  and not held in registers,  or
*        are not affected by optimizers.
*
*        DLAMC1 returns the parameters  LBETA, LT, LRND and LIEEE1.
*
         CALL DLAMC1( LBETA, LT, LRND, LIEEE1 )
*
*        Start to find EPS.
*
         B = LBETA
         A = B**( -LT )
         LEPS = A
*
*        Try some tricks to see whether or not this is the correct  EPS.
*
         B = TWO / 3
         HALF = ONE / 2
         SIXTH = DLAMC3( B, -HALF )
         THIRD = DLAMC3( SIXTH, SIXTH )
         B = DLAMC3( THIRD, -HALF )
         B = DLAMC3( B, SIXTH )
         B = ABS( B )
         IF( B.LT.LEPS )
     $      B = LEPS
*
         LEPS = 1
*
*+       WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
   10    CONTINUE
         IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
            LEPS = B
            C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
            C = DLAMC3( HALF, -C )
            B = DLAMC3( HALF, C )
            C = DLAMC3( HALF, -B )
            B = DLAMC3( HALF, C )
            GO TO 10
         END IF
*+       END WHILE
*
         IF( A.LT.LEPS )
     $      LEPS = A
*
*        Computation of EPS complete.
*
*        Now find  EMIN.  Let A = + or - 1, and + or - (1 + BASE**(-3)).
*        Keep dividing  A by BETA until (gradual) underflow occurs. This
*        is detected when we cannot recover the previous A.
*
         RBASE = ONE / LBETA
         SMALL = ONE
         DO 20 I = 1, 3
            SMALL = DLAMC3( SMALL*RBASE, ZERO )
   20    CONTINUE
         A = DLAMC3( ONE, SMALL )
         CALL DLAMC4( NGPMIN, ONE, LBETA )
         CALL DLAMC4( NGNMIN, -ONE, LBETA )
         CALL DLAMC4( GPMIN, A, LBETA )
         CALL DLAMC4( GNMIN, -A, LBETA )
         IEEE = .FALSE.
*
         IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
            IF( NGPMIN.EQ.GPMIN ) THEN
               LEMIN = NGPMIN
*            ( Non twos-complement machines, no gradual underflow;
*              e.g.,  VAX )
            ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
               LEMIN = NGPMIN - 1 + LT
               IEEE = .TRUE.
*            ( Non twos-complement machines, with gradual underflow;
*              e.g., IEEE standard followers )
            ELSE
               LEMIN = MIN( NGPMIN, GPMIN )
*            ( A guess; no known machine )
               IWARN = .TRUE.
            END IF
*
         ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
            IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
               LEMIN = MAX( NGPMIN, NGNMIN )
*            ( Twos-complement machines, no gradual underflow;
*              e.g., CYBER 205 )
            ELSE
               LEMIN = MIN( NGPMIN, NGNMIN )
*            ( A guess; no known machine )
               IWARN = .TRUE.
            END IF
*
         ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND.
     $            ( GPMIN.EQ.GNMIN ) ) THEN
            IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
               LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
*            ( Twos-complement machines with gradual underflow;
*              no known machine )
            ELSE
               LEMIN = MIN( NGPMIN, NGNMIN )
*            ( A guess; no known machine )
               IWARN = .TRUE.
            END IF
*
         ELSE
            LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
*         ( A guess; no known machine )
            IWARN = .TRUE.
         END IF
***
* Comment out this if block if EMIN is ok
         IF( IWARN ) THEN
            FIRST = .TRUE.
            WRITE( 6, FMT = 9999 )LEMIN
         END IF
***
*
*        Assume IEEE arithmetic if we found denormalised  numbers above,
*        or if arithmetic seems to round in the  IEEE style,  determined
*        in routine DLAMC1. A true IEEE machine should have both  things
*        true; however, faulty machines may have one or the other.
*
         IEEE = IEEE .OR. LIEEE1
*
*        Compute  RMIN by successive division by  BETA. We could compute
*        RMIN as BASE**( EMIN - 1 ),  but some machines underflow during
*        this computation.
*
         LRMIN = 1
         DO 30 I = 1, 1 - LEMIN
            LRMIN = DLAMC3( LRMIN*RBASE, ZERO )
   30    CONTINUE
*
*        Finally, call DLAMC5 to compute EMAX and RMAX.
*
         CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
      END IF
*
      BETA = LBETA
      T = LT
      RND = LRND
      EPS = LEPS
      EMIN = LEMIN
      RMIN = LRMIN
      EMAX = LEMAX
      RMAX = LRMAX
*
      RETURN
*
 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
     $      '  EMIN = ', I8, /
     $      ' If, after inspection, the value EMIN looks',
     $      ' acceptable please comment out ',
     $      / ' the IF block as marked within the code of routine',
     $      ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / )
*
*     End of DLAMC2
*
      END
*
************************************************************************
*
      DOUBLE PRECISION FUNCTION DLAMC3( A, B )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION   A, B
*     ..
*
*  Purpose
*  =======
*
*  DLAMC3  is intended to force  A  and  B  to be stored prior to doing
*  the addition of  A  and  B ,  for use in situations where optimizers
*  might hold one of these in a register.
*
*  Arguments
*  =========
*
*  A, B    (input) DOUBLE PRECISION
*          The values A and B.
*
* =====================================================================
*
*     .. Executable Statements ..
*
      DLAMC3 = A + B
*
      RETURN
*
*     End of DLAMC3
*
      END
*
************************************************************************
*
      SUBROUTINE DLAMC4( EMIN, START, BASE )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. Scalar Arguments ..
      INTEGER            BASE, EMIN
      DOUBLE PRECISION   START
*     ..
*
*  Purpose
*  =======
*
*  DLAMC4 is a service routine for DLAMC2.
*
*  Arguments
*  =========
*
*  EMIN    (output) EMIN
*          The minimum exponent before (gradual) underflow, computed by
*          setting A = START and dividing by BASE until the previous A
*          can not be recovered.
*
*  START   (input) DOUBLE PRECISION
*          The starting point for determining EMIN.
*
*  BASE    (input) INTEGER
*          The base of the machine.
*
* =====================================================================
*
*     .. Local Scalars ..
      INTEGER            I
      DOUBLE PRECISION   A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMC3
      EXTERNAL           DLAMC3
*     ..
*     .. Executable Statements ..
*
      A = START
      ONE = 1
      RBASE = ONE / BASE
      ZERO = 0
      EMIN = 1
      B1 = DLAMC3( A*RBASE, ZERO )
      C1 = A
      C2 = A
      D1 = A
      D2 = A
*+    WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
*    $       ( D1.EQ.A ).AND.( D2.EQ.A )      )LOOP
   10 CONTINUE
      IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND.
     $    ( D2.EQ.A ) ) THEN
         EMIN = EMIN - 1
         A = B1
         B1 = DLAMC3( A / BASE, ZERO )
         C1 = DLAMC3( B1*BASE, ZERO )
         D1 = ZERO
         DO 20 I = 1, BASE
            D1 = D1 + B1
   20    CONTINUE
         B2 = DLAMC3( A*RBASE, ZERO )
         C2 = DLAMC3( B2 / RBASE, ZERO )
         D2 = ZERO
         DO 30 I = 1, BASE
            D2 = D2 + B2
   30    CONTINUE
         GO TO 10
      END IF
*+    END WHILE
*
      RETURN
*
*     End of DLAMC4
*
      END
*
************************************************************************
*
      SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. Scalar Arguments ..
      LOGICAL            IEEE
      INTEGER            BETA, EMAX, EMIN, P
      DOUBLE PRECISION   RMAX
*     ..
*
*  Purpose
*  =======
*
*  DLAMC5 attempts to compute RMAX, the largest machine floating-point
*  number, without overflow.  It assumes that EMAX + abs(EMIN) sum
*  approximately to a power of 2.  It will fail on machines where this
*  assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
*  EMAX = 28718).  It will also fail if the value supplied for EMIN is
*  too large (i.e. too close to zero), probably with overflow.
*
*  Arguments
*  =========
*
*  BETA    (input) INTEGER
*          The base of floating-point arithmetic.
*
*  P       (input) INTEGER
*          The number of base BETA digits in the mantissa of a
*          floating-point value.
*
*  EMIN    (input) INTEGER
*          The minimum exponent before (gradual) underflow.
*
*  IEEE    (input) LOGICAL
*          A logical flag specifying whether or not the arithmetic
*          system is thought to comply with the IEEE standard.
*
*  EMAX    (output) INTEGER
*          The largest exponent before overflow
*
*  RMAX    (output) DOUBLE PRECISION
*          The largest machine floating-point number.
*
* =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
*     ..
*     .. Local Scalars ..
      INTEGER            EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
      DOUBLE PRECISION   OLDY, RECBAS, Y, Z
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMC3
      EXTERNAL           DLAMC3
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MOD
*     ..
*     .. Executable Statements ..
*
*     First compute LEXP and UEXP, two powers of 2 that bound
*     abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
*     approximately to the bound that is closest to abs(EMIN).
*     (EMAX is the exponent of the required number RMAX).
*
      LEXP = 1
      EXBITS = 1
   10 CONTINUE
      TRY = LEXP*2
      IF( TRY.LE.( -EMIN ) ) THEN
         LEXP = TRY
         EXBITS = EXBITS + 1
         GO TO 10
      END IF
      IF( LEXP.EQ.-EMIN ) THEN
         UEXP = LEXP
      ELSE
         UEXP = TRY
         EXBITS = EXBITS + 1
      END IF
*
*     Now -LEXP is less than or equal to EMIN, and -UEXP is greater
*     than or equal to EMIN. EXBITS is the number of bits needed to
*     store the exponent.
*
      IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
         EXPSUM = 2*LEXP
      ELSE
         EXPSUM = 2*UEXP
      END IF
*
*     EXPSUM is the exponent range, approximately equal to
*     EMAX - EMIN + 1 .
*
      EMAX = EXPSUM + EMIN - 1
      NBITS = 1 + EXBITS + P
*
*     NBITS is the total number of bits needed to store a
*     floating-point number.
*
      IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
*
*        Either there are an odd number of bits used to store a
*        floating-point number, which is unlikely, or some bits are
*        not used in the representation of numbers, which is possible,
*        (e.g. Cray machines) or the mantissa has an implicit bit,
*        (e.g. IEEE machines, Dec Vax machines), which is perhaps the
*        most likely. We have to assume the last alternative.
*        If this is true, then we need to reduce EMAX by one because
*        there must be some way of representing zero in an implicit-bit
*        system. On machines like Cray, we are reducing EMAX by one
*        unnecessarily.
*
         EMAX = EMAX - 1
      END IF
*
      IF( IEEE ) THEN
*
*        Assume we are on an IEEE machine which reserves one exponent
*        for infinity and NaN.
*
         EMAX = EMAX - 1
      END IF
*
*     Now create RMAX, the largest machine number, which should
*     be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
*
*     First compute 1.0 - BETA**(-P), being careful that the
*     result is less than 1.0 .
*
      RECBAS = ONE / BETA
      Z = BETA - ONE
      Y = ZERO
      DO 20 I = 1, P
         Z = Z*RECBAS
         IF( Y.LT.ONE )
     $      OLDY = Y
         Y = DLAMC3( Y, Z )
   20 CONTINUE
      IF( Y.GE.ONE )
     $   Y = OLDY
*
*     Now multiply by BETA**EMAX to get RMAX.
*
      DO 30 I = 1, EMAX
         Y = DLAMC3( Y*BETA, ZERO )
   30 CONTINUE
*
      RMAX = Y
      RETURN
*
*     End of DLAMC5
*
      END
      LOGICAL          FUNCTION LSAME( CA, CB )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     September 30, 1994
*
*     .. Scalar Arguments ..
      CHARACTER          CA, CB
*     ..
*
*  Purpose
*  =======
*
*  LSAME returns .TRUE. if CA is the same letter as CB regardless of
*  case.
*
*  Arguments
*  =========
*
*  CA      (input) CHARACTER*1
*  CB      (input) CHARACTER*1
*          CA and CB specify the single characters to be compared.
*
* =====================================================================
*
*     .. Intrinsic Functions ..
      INTRINSIC          ICHAR
*     ..
*     .. Local Scalars ..
      INTEGER            INTA, INTB, ZCODE
*     ..
*     .. Executable Statements ..
*
*     Test if the characters are equal
*
      LSAME = CA.EQ.CB
      IF( LSAME )
     $   RETURN
*
*     Now test for equivalence if both characters are alphabetic.
*
      ZCODE = ICHAR( 'Z' )
*
*     Use 'Z' rather than 'A' so that ASCII can be detected on Prime
*     machines, on which ICHAR returns a value with bit 8 set.
*     ICHAR('A') on Prime machines returns 193 which is the same as
*     ICHAR('A') on an EBCDIC machine.
*
      INTA = ICHAR( CA )
      INTB = ICHAR( CB )
*
      IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN
*
*        ASCII is assumed - ZCODE is the ASCII code of either lower or
*        upper case 'Z'.
*
         IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32
         IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32
*
      ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN
*
*        EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
*        upper case 'Z'.
*
         IF( INTA.GE.129 .AND. INTA.LE.137 .OR.
     $       INTA.GE.145 .AND. INTA.LE.153 .OR.
     $       INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64
         IF( INTB.GE.129 .AND. INTB.LE.137 .OR.
     $       INTB.GE.145 .AND. INTB.LE.153 .OR.
     $       INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64
*
      ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN
*
*        ASCII is assumed, on Prime machines - ZCODE is the ASCII code
*        plus 128 of either lower or upper case 'Z'.
*
         IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32
         IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32
      END IF
      LSAME = INTA.EQ.INTB
*
*     RETURN
*
*     End of LSAME
*
      END
C----------------------------------------------------------------------------
